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This review is focused on the borderhne region of theoretical physics and mathemat- 
ics. First, we describe numerical methods for the acceleration of the convergence of 
series. These provide a useful toolbox for theoretical physics which has hitherto not 
received the attention it actually deserves. The unifying concept for convergence 
acceleration methods is that in many cases, one can reach much faster convergence 
than by adding a particular series term by term. In some cases, it is even possible 
to use a divergent input series, together with a suitable sequence transformation, 
for the construction of numerical methods that can be applied to the calculation of 
special functions. This review both aims to provide some practical guidance as well 
as a groundwork for the study of specialized literature. As a second topic, we re- 
view some recent developments in the field of Borel resummation, which is generally 
recognized as one of the most versatile methods for the summation of factorially 
divergent (perturbation) series. Here, the focus is on algorithms which make opti- 
mal use of all information contained in a finite set of perturbative coefBcients. The 
unifying concept for the various aspects of the Borel method investigated here is 
given by the singularities of the Borel transform, which introduce ambiguities from 
a mathematical point of view and lead to different possible physical interpretations. 
The two most important cases are: (i) the residues at the singularities correspond 
to the decay width of a resonance, and (ii) the presence of the singularities indi- 
cates the existence of nonperturbative contributions which cannot be accounted for 
on the basis of a Borel resummation and require generalizations toward resurgent 
expansions. Both of these cases are illustrated by examples. 
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1 Background and Orientation 



The purpose of this section is that of being a starting point and a reference for the reader. Its 
contents are organized as follows: 

• In Sec. II. H a general introduction to convergence acceleration and resummation is provided. 

• In Sec. \1.2\ some basic terminology is introduced. This terminology is a basis for the general 
familiarity with the notions that govern the construction of convergence acceleration methods. 

The concepts contained in this section constitute to a large extent a prerequisite for the material 
presented in the following sections. 

1.1 Introduction to Convergence Acceleration and Resummation 

The basic idea of convergence acceleration is to extract hidden information contained in partial sums 
of a specific slowly convergent series, and to use that information in order to make, loosely speaking, 
an educated guess about higher-order partial sums which eventually converge to some limit. In 
many cases, the "educated guesses" can lead to spectacular numerical results which represent a 
drastic improvement over a term-by-term summation of the series, even if the series is formally 
convergent. 

Convergence acceleration methods have a long history. E.g., as mentioned of p. ix of Ref. [1], such 
methods were already used by the Chinese mathematicians Liu Hui (A.D. 263) and Zhu Chongzhi 
(429 — 500) for obtaining better approximations to vr. It is thus clear that, at least in rudimentary 
form, convergence acceleration methods are much older than calculus. Another reference to early 
successful attempts can be found on pp. 90 — 91 of a book by Brezinski [2] where it is mentioned 
that convergence acceleration methods have already been used in 1654 by Huygens and in 1674 
by Seki Kowa, the probably most famous Japanese mathematician of his time. Both Huygens and 
Seki Kowa were trying to obtain better approximations to vr. Huygens used a linear extrapolation 
scheme which is a special case of what we now call Richardson extrapolation [3], and Seki Kowa 
used the so-called process, which is usually attributed to Aitken [4]. 

Sequence transformations are principal numerical tools to overcome convergence problems. In this 
approach, a slowly convergent or divergent sequence {sn}^=Q, whose elements may for instance be 
the partial sums s„ = X^^^q of an infinite series, is converted into a new sequence {s'„}^q 
with hopefully better numerical properties. Before the invention of electronic computers, mainly 
linear sequence transformations were used. These matrix transformations compute the elements 
of the transformed sequence {s'„}^q as weighted averages of the elements of the input sequence 
{sn}'^=Q according to the general scheme s'„ = X]fc=o l^nkSk- The theoretical properties of such 
matrix transformations are very well understood [5-13]. Their main appeal lies in the fact that 
for the weights some necessary and sufficient conditions could be formulated which guarantee 
that the application of such a matrix transformation to a convergent sequence {s„}^q yields a 
transformed sequence {s^}^q converging to the same limit s = Soo- Theoretically, this regularity 
is extremely desirable, but from a practical point of view, it is a disadvantage. This may sound 
paradoxical. However, according to a remark on p. X of Ref. [10], the size of the domain of regularity 
of a transformation and its efficiency seem to be inversely related. Accordingly, regular linear 
transformations are in general at most moderately powerful, and the popularity of most linear 
transformations has declined considerably in recent years. 

Modern nonlinear sequence transformations as for instance the e algorithm [14] and the p algo- 
rithm [15] or the i9 algorithm [16] have largely complementary properties: They are nonregular. 
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which means that that the convergence of the transformed sequence is not guaranteed, let alone to 
the correct limit. In addition, the theoretical properties of nonlinear transformations are far from 
being completely understood. Nevertheless, they often accomplish spectacular results when ap- 
plied to appropriate problems. Consequently, nonlinear transformations now dominate both math- 
ematical research as well as practical applications, as documented in a large number of recent 
books [1,10,17-29] and review articles [30-33] on this topic. 

We here attempt to present these transformations in both a special context (with a focus on 
applications) as well as in the general context given by an underlying principle of construction 
exemplified by the E algorithm. One of the recurrent themes in the construction of sequence 
transformations is the separation of the series into a partial sum s„ = X^^^g ^^'^ ^ remainder 
Tn = — Yl'h=n+i SO that s„, = s + r„. The remainder term r„ is to be approximated by remainder 
estimates and eliminated from the input data s„ by suitable transformations. 

We restrict the discussion to the acceleration of sequences. Quite naturally, convergence accelera- 
tion methods can be proposed and applied to the pointwise evaluation of partial-wave expansions 
which occur ubiquitously in physics (see, e.g., Ref. [34,35]), and to Fourier expansions and other 
orthogonal series [36]. Here, the sequences are typically composed of the partial sums of the re- 
spective expansions. For applications to vector or matrix problems, the reader is referred to the 
books [25,37] and numerous articles, e.g. [38,39], and for applications to numerical quadrature, in 
particular with respect to oscillatory functions to be integrated over semiinfinite intervals, we refer 
to the books [28,40] and several articles which describe corresponding applications, e.g. [41]. 

Forming a bridge between convergence acceleration and resummation, one may observe that some 
acceleration algorithms may also be used for the resummation of (mildly) divergent series, but are 
in general less versatile than the Borel method. The resummed value of a divergent series, which is 
sometimes called antilimit, can be attributed, according to Euler, to the particular finite mathemat- 
ical entity whose expansion gives rise to the original series. L. Euler wrote in a letter to Goldbach 
(1745): 'Summa cuiusque seriei est valor expressionis illius finitae, ex cuius evolutione ilia series 
oritur, " which may be translated from the Latin language as: "The sum of any given series [whether 
convergent or divergent] is the value of the specific finite expression whose expansion gave rise to 
that same series." One can mention a number of accepted methods for series (re-)summation: e.g., 
there are averaging processes for oscillatory partial sums like the Cesaro method [6] or sophisticated 
mathematical constructions like zeta function regularization [42] , which depend on the availability 
of closed-form analytic expressions for all terms that need to be summed, and permit the variation 
of a free parameter (the "argument" of the generalized zeta function) from the physically relevant 
to a mathematically favourable domain, where the defining series becomes convergent. Later, in 
order to allow for the determination of the physical property of interest, one performs an analytic 
continuation of the parameter back to the physical domain. In the same context, one can mention 
the Mellin-Barnes transform [43]. 

Divergent series result from any eigenvalue perturbation theory, unless the perturbation is Kato- 
bounded with respect to the unperturbed Hamiltonian (see Ref. [44]). In particular, the application 
of simple, "naive" Rayleigh-Schrodinger perturbation theory to a harmonic oscillator with an x'^ 
perturbation (where x is the position operator) already leads to a divergent series in higher orders 
of perturbation theory. While the degree of acceptance of divergent series within the mathematical 
scientific community underwent some undeserved oscillations during the last two centuries, we are 
today in the fortunate situation that their properties are routinely investigated and harvested in 
the context of physical problems that originate from perturbation theory in quantum mechanics 
and field theory. 

One of the most versatile and most relevant resummation algorithms, especially those originating 
from perturbation theory, is the Borel method [45-47]. The reason is that many physical diver- 
gences are factorial, and for these, the Borel summation is the most powerful one. Already in 1952 
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Dyson [48] had argued that perturbation expansions in quantum electrodynamics should diverge, 
and it later turned out [49-54] that the divergence should be factorial. Also, Bender and Wu [55-57] 
demonstrated that Dyson's speculation was not a mere mathematical sophistication: They proved 
that the coefficients of the standard Rayleigh-Schrodinger perturbation expansions for the energy 
eigenvalues of the anharmonic oscillators with anharmonicities x^'" with m> 2 grow at least facto- 
rially with the index n as (m n) ! . A new subfield of theoretical physics emerged which is nowadays 
called large-order perturbation theory (see for example [58-60] and references therein). Following 
these works, many other quantum mechanical systems described by perturbation expansions were 
investigated, and in the overwhelming majority of all cases very similar divergence problems were 
observed. Moreover, nowadays the factorial divergence of perturbation expansions in quantum field 
theory and in related disciplines now almost has the status of a law of nature (see for example the 
discussion in [61]). 

The Borel resummation process replaces the divergent series, after analytic continuation of the 
Borel transform, by a Laplace-type integral, which is sometimes called the Laplace-Borel integral 
for obvious reasons. A famous ambiguity results in the case of singularities of the transform along 
the integration contour. In principle, one may imagine three possible interpretations for the am- 
biguity introduced by the singularities, and it is interesting to remark that all of them might be 
realized in nature. They are marked here as (i), (ii) and (iii) and will be briefly discussed now. 
(i) In some cases (e.g., the cubic anharmonic oscillator and the Stark effect), the resummation is 
achieved by using a complex integration contour, which goes around the singularity on a half circle 
and leads to an imaginary part for the (resonance) energy eigenvalues of the Hamiltonians under 
investigation. Of course, the sign of the imaginary part depends on the clockwise or anticlockwise 
circulation around the pole and cannot be determined on the basis of the purely real perturbation 
series alone: additional physics input is needed and provided by the very nature of a resonance. 
These cases find a natural mathematical interpretation via the concept of distributional Borel 
summability. (ii) There are situations in which the singularities of the Borel transform indicate 
the presence of further contributions that cannot be obtained on the basis of perturbation theory 
alone, even if perturbation theory is combined with the most sophisticated resummation methods. 
Rather, one has to generalize perturbation theory to a resurgent expansion [62-65]. In some cases, 
there are subtle cancellations of the imaginary parts created by the singularities of the Borel trans- 
form of the perturbation series against explicit imaginary parts of higher-order instanton effects 
which are characterized by nonanalytic factors, leading to purely real eigenenergies for, e.g., quan- 
tum mechanical potentials of the double-well type. The instanton effects, in turn, are decorated 
again by divergent nonalternating perturbation series, which give to yet further imaginary parts, 
but these are again compensated by instanton effects of even higher order. The generalized per- 
turbative expansion which entails all of these effects is a so-called resurgent expansion, (iii) One 
may easily construct mathematical model problems where the correct resummation is achieved by 
a principal- value prescription for the evaluation of the Borel integral. It has been conjectured that 
this prescription might be the correct one for a number of physical problems as well. Summarizing, 
we may state that the ambiguity of the Borel resummation process definitely exists. This ambiguity 
provides for an interpretation of the Borel transform even in cases where the complete physical 
answer cannot be obtained on the basis of perturbation theory alone. 

The resummation methods allow us to determine eigenenergies and other physical properties on 
the basis of a finite number of perturbative coefficients alone. One can now even go one step further 
and observe that, on the one hand, the large-order behaviour of perturbation theory is connected to 
the underlying perturbative potential or interaction Lagrangian (in the case of a field theory), and 
yet, on the other hand, the potential also determines the large-coupling expansion for the potential 
or the field theory under investigation. So, one can even try to guess the large-coupling expansion 
from a finite number of terms of the weak-coupling perturbative expansion, because these already 
give an impression of the large-order behaviour. This paradoxical endeavour has been pursued in 
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recent years, in a series of papers, by Suslov [66-69]. 

Let us now attempt to form a bridge between convergence acceleration and resummation via se- 
quence transformations, which can be used for the resummation of (mildly) divergent series in 
addition to their use as convergence accelerators. We recall that the convergence acceleration of 
the sequence {s„}^q is achieved by eliminating the remainders r„ from a necessarily finite set 
of input data {sn}^=Q, in order to obtain a better estimate s'^ to s than the obvious estimate 
which would be provided by Sm- Clearly, the "elimination of remainders" is a concept which can 
be transported to the case of a divergent series. The only difference is that for a convergent series, 
^ as n — > oo, whereas for a divergent series, lim^^oo does not exist. In some cases, the 
sequence transformations provide us with numerically favourable, highly efficient algorithms for 
the computation of the antilimit of the divergent input series, so that they constitute competitive 
numerical algorithms for the actual computation of the antilimit, which could otherwise only ob- 
tained by direct evaluation of the expression which gave rise to the divergent series in the first 
place, or by the Borel method in connection with Pade approximants for the analytic continuation 
of the Borel transform. The latter method, although versatile, is numerically less efficient because 
it entails, at least, the numerical evaluation of the Laplace-Borel integral. 



1.2 Mathematical Basis 

1.2.1 Special Mathematical Symbols 

In order to fix ideas, it is perhaps useful to recall a few mathematical symbols. In particular, 

N = {1,2,3,...} (1) 

denotes the set of positive integers. No = N U {0} stands for the set {0, 1, 2, . . .} of nonnegative 
integers, and 

Z = {0,±1,±2,...} (2) 

is the set of positive and negative integers. Moreover, M and C denote the sets of real and complex 
numbers, respectively. 

The symbol is the integral part of x G M, i.e., the largest integer m £ Z satisfying the inequality 
m < X. 

The nonstandard notation 

{Sn}'^=0 (3) 

is used to denote a sequence of elements Sn with n £ No- In this work, we usually assume that 
every element Sn represents the partial sum of an infinite series: 

n 

Indeed, such a definition implies that the sequence elements s_i, s_2) >5_3, . . . with negative indices 
are empty sums and thus equal to zero. 

As just indicated, empty sums are always interpreted as zero: 

n 

a/c = , if m > n . (5) 

k=in 

Similarly, an empty product is always interpreted as one: 

n 

JJ afc = 1 , if m > n . (6) 

k=m 



?n = ^ak. (4) 
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Let f{z) and g{z) be two functions defined on some domain D in the complex plane and let zq be 
a limit point of D, possibly the point at infinity. Then, 

f{z) = 0{g{z)) , z^zo, (7) 

means that there is a positive constant A and a neighborhood U of zq such that 

1/(^)1 < A\g{z)\ (8) 

for all z £ U D D. If g{z) does not vanish on U f) D this simply means that f{z)/g{z) is bounded 
on U D D. Similarly, 

f{z) = o{g{z)) , z^zo, (9) 
means that for any positive number e G ffi there exists a neighborhood U of zq such that 

1/(^)1 < (10) 

for all z € U n D. If g(z) does not vanish on U f] D this simply means that f{z)/g{z) approaches 
zero as z ^ Zq. 

1.2.2 Finite Differences 

So-called finite difference operators [70, 71] are essential in the construction of sequence transfor- 
mations. Let us consider a function /(n) defined on the set No of nonnegative integers. Then, the 
forward difference operator A/(n) is defined by 

A/(n) = /(n+l)-/(n), n G Nq . (11) 

Higher powers of this operator can be defined recursively: 

^^f{n) = ^[^^-^f{n)], fcGN, (12) 

with A''/(n) = /(n). Apart from the forward difference one may also introduce the shift operator: 

Ef{n) = /(n + 1), (13) 

whose higher powers can again be defined recursively: 

E^fin) = E[E^-^f{n)] = f{n + k), A:eN, (14) 

with E^f{n) = f{n). As seen from definitions (jlip and ()13p . the forward difference A and shift E 
operators are connected by the relationship 

A = E-l, (15) 

where 1 is the identity operator. This relationship can be combined with the binomial theorem to 
give 

AV(n) = (-1)^ fl (-1)' C) + ■?■)' ^ e No . (16) 

In the following, we always tacitly assume that in the case of quantities, which depend on several 
indices like the difference operator A and the shift operator E only act only upon n and not 
on any other indices: e.g., we have A^^""^ = ^^""''^^ — .4.[,"\ 

Furthermore, we should clarify the notation for sequences with /(n) = Sn- The quantity As„ is to 
be understood as As„ = s„+i — Sn = fln+i, and we have A^s„ = Sn+2 — 2s„+i + Sn = Aa„_|_i = 
an+2 — CLn+i and so on [we assume the structure given in Eq. (jH) for the Sn]- 
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1.2.3 Asymptotic Sequences and Asymptotic Expansions 

A finite or infinite sequence of functions {^niz)}^^Q, whicli are defined on some domain D of 
complex numbers on wliich all ^niz) are nonzero except possibly at zq, is called an asymptotic 
sequence as z ^ zq if for all n G No, 

$n+i(z) = o(^>„(z)), z^zo- (17) 

The formal series 

oo 

/(z)~^a„cl>„(z), (18) 

n=0 

which need not be convergent, is called an asymptotic expansion of f{z) with respect to the 
asymptotic sequence {^n{z)}^^Q in the sense of Poincare [72] if, for every m G No, 

m 

/(z) - ^a„$„(z) = o{<^„^{z)), Z^Zq. (19) 

n=0 

If such a Poincare-type asymptotic expansion exists, it is unique, and its coefficients a„ can be 
computed as follows, 

m—l 

am = ^lim I [f{z) - an^niz)] I «>m(z)} , m G No . (20) 

n=0 

The particularly important case of an asymptotic expansion with respect to the asymptotic se- 
quence {2;~"'}^o about = oo has the form of an asymptotic power series in 1/z, 

oo 

f[z) ~ "fc/^''- z^<x,. (21) 

fc=0 

In the current article, we will often use the notation 

oo 

/(z) ~ J^Ofc/, ^^0, (22) 

fc=0 

for asymptotic expansions or perturbative expansions of a physical quantity /(z) with respect to 
the asymptotic sequence {^^"'l^o ™ sense of Poincare. 

1.2.4 Remainder Estimates 

Let us assume that {s„}^o either converges to some limit s, or, if it diverges, can be summed by 
an appropriate summation method to give s. Following [73], this generalized limit s is frequently 
called antilimit in the case of divergence. In either case, there is a partition of a sequence element 
Sn into the limit or antilimit s and the remainder r„ according to 

Sn = S ^ Tn (23) 

for all n G Nq. If s„ is the partial sum 

n 

Sn = Yak (24) 

fe=0 
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of an infinite series, then the remainder obviously satisfies 

oo 

r„ = - ^ ak- (25) 

fc=n+l 

In particular, this implies that the remainders r„ are all negative if the series terms an are all 
positive. 

The dependence of the magnitude of the remainder on the index n is a natural measure for the 
convergence of the sequences or series. Often, it is of considerable interest to analyze the asymptotics 
of the sequence of remainders {r„}^Q as n ^ oo. Let us assume that we have found an asymptotic 
sequence {<^fc(n.)}|^Q as n ^ oo with (po{n) = 1 for the function rn/tOn, which is equivalent to 
identifying some dominant part (leading term) cvn of with respect to the asymptotic sequence 
{(^fc(n)}^Q. We call cOn the remainder estimate. In particular, this means that the remainder term 
r„, can be expanded asymptotically according to 

oo 

rn/^n ~ X] '^'^ '^^^^'^ ' OO. (26) 

fc=0 

Sequences of remainder estimates {ujn}^=Q will be of considerable importance in this report. The 
reason is that it is often possible to obtain at least some structural information about the behavior 
of the dominant term of a remainder as n ^ co. It will become clear later that convergence accel- 
eration or summation methods, which explicitly utilize the information contained in the remainder 
estimates {u;„}^0) frequently lead to efficient algorithms. 



1.2.5 Types of Convergence 

The ideas of "linear" and "logarithmic" convergence turn out to be very useful concepts when 
analyzing the types of convergence of series, and are used very often in the literature on convergence 
acceleration. They are recalled in the current section. 

We observe that the behaviour of the remainder (j25p as a function of n is a natural measure for 
the rate of convergence of the series. In this subsection, therefore, we discuss how different types 
of behaviour of the r„ at n — > oo may be used to classify the types of convergence of a particular 
sequence. In order to perform such a classification, we write 

lim £n±i^ = lim !ii±i = p, (27) 

n— »oo Sn — S n^oo r^j 

where, of course, we assume that the sequence {snlj^o converges to some limit s. If, in this equation, 
< |/o| < 1, we say that the sequence {s^j^o converges linearly, if p = 1 holds, we say that this 
sequence converges logarithmically, and if /? = holds, we say that it converges hyper linearly. 
Moreover, \p\ > 1 implies that the sequence {s„}^o diverges. 

In Ref . [74] , it has been shown that if the terms of a convergent series are all real and have the 
same sign, then for s„ = X]fc=o^'=' ^™ establish a connection between the remainder r„ and 
the series term Qjifi clS follows, 

lim !j1±1 = Ihn ^ = p. (28) 

n— >oo Tn n^oo CL^ 

Effectively, this means that we can replace in the fraction r„+i/r„ each of the r„ by As^, 

lim "^"^ = lim = p . (29) 

n— >oo Tn n— >oo As»j 
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Apart from the particular case of the real and nonalternating a^, Eq. (j29p can be also applied to 
arbitrary hyperlinearly and linearly convergent real sequences [10]. 

We are ready now to illustrate the application of Eq. (j27p for the classification of the (type of) 
convergence of some particular series. To this end, let us start from the Riemann zeta function, 
defined by a Dirichlet series as given in Ref. [75], 

((z) = Y — . (30) 

n=0 ^ ^ 

We would like to determine the type of convergence of this series for a particular case of z = 2. 
For this case, the reminder (1251) reads: 



oo ^ 

pi) 

=n+l ^ 



By applying the Euler-Maclaurin formula [see, e.g., Eqs. (2.01) — (2.02) on p. 285 of Ref. [76]] to 
the right side of this expression, we find the leading and the first sub-leading asymptotics of the 
reminder r„, 

1 1 ^. 



n + 2 2(n + 2)2 ^ (n + 2)2i+i 

Here, B2j is a Bernoulli number [p. 281 of Ref. [76]]. Together with Eq. (127p . this immediately 
implies that the series pop at z = 2 converges logarithmically, i.e. p = 1 for n ^ oo. 

A more elaborate example concerns the convergence properties of the alternating series with partial 
sums: 

" (_l)fc 7r2 
k=0 ^ ' 

This series is similar to the Dirichlet series for the zeta function C(2), but is strictly alternating 
in sign. The fact that it converges to -k^ jVl is shown, for example, in Sec. 2.2 of Ref. [77]. Via a 
somewhat rather lengthy calculation, we find for ?i ^ cxo: 

^ (-1)'= (-1)" / I 3 6 9 ^, 

By inserting this reminder into Eq. (j27p . we obtain that p = —\. This means that while \p\ = 1, we 
have p ^ ^ and so the series is not logarithmically convergent. In this contribution, we will refer 
to convergent series with remainder terms r„_|_i/r„ — > —1 for n — > oo (i.e., series with p = —1) as 
alternating marginally convergent series. 

So far, we have discussed the application of the remainder estimates for the classification of the 
convergence of a particular single sequence. The estimates r„, however, can be also be applied 
to compare the rate of convergence of two different series. To this end, let us assume that two 
sequences {sn}^=Q and {s'n\^=Q both converge to the same limit s. We shall say that the sequence 
{s^jJ^Q converges more rapidly than {s„}^q if 

lim = lim ^ = 0. (35) 

n^oo Sn — S n— >oo 

A further natural question is whether one can replace in Eq. (j35]) r„ by As.„ and r'^ by As'„, in 
analogy to Eq. (P^|) : 

lim f!i±l^ = lim ^ = 0. (36) 



n— >oo 



Sn+l - Sn ASn 
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Note that in typical cases, the quantity s'^_^i — s'^ which may result from the application of a 
sequence transformation to {sn}'^=Q cannot simply be replaced by some known term a„. 

It seems (see Ref. [31]) that the equivalence of ([35]) and ([36]) cannot be established universally 
without making some explicit assumptions about how fast the sequences {s„}J^q and {s^}^q 
approach their common limit s. In particular, if {s„}^q converges linearly, then the transformed 
series {s^}^g can converge more rapidly only if it converges at least linearly or even faster. We 
can always write 

^ <+i - 4 ^ <-g - g) / (4 - g)] - 1 ^ - 1 ^^^^ 

ASn Sn+l - Sn Sn - S [{Sn+1 - s) / (s„ - s)] - 1 r„+i/r„ - 1 ' 

Let us assume now that both {s',j}^q as well as {sn}.^o fulfill Eq. ([29|) with limits p' and p. We 
can then write 

lim ^ = hm ^ . (38) 

n— >oo As„ n^oo r,^ p — 1 

The equivalence of (j35p and (j36p immediately follows under the assumptions made, which include 
the linear convergence of the input series (this ensures p ^ 1). If, however, {s„}^q converges 
logarithmically, the denominator of the second term on the right-hand side of Eq. (j38p approaches 
zero as n — > oo. In this case, some additional assumptions about the rate of convergence of {s^}^q 
to s have to be made. 

Finally, we should perhaps remark that the analogous terminology of "exponential" and "algebraic" 
convergence (instead of "linear" and "logarithmic" convergence) is often used in both the physical 
as well as the mathematical literature, in particular in situations where the "convergent entities" 
are different from partial sums of a series. An example would be the rate of convergence of iterative 
solutions to dynamical problems. 



1.2.6 Stieltjes Functions and Stieltjes Series 



Stieltjes series are, roughly speaking, paradigmatic examples of asymptotic series for which it is 
possible to find remainder estimates. These estimates characterize the remainder after the calcula- 
tion of a partial sum and can be used for the construction of sequence transformations which profit 
from the additional information contained in the remainder estimates (see, for example. Section II 
in [78] and references therein). 

Specifically, a function F{z) with z G C is called a Stieltjes function if it can be expressed as a 
Stieltjes integral according to 

= |arg(z)|<vr. (39) 

Here, ^{t) is a bounded, nondecreasing function taking infinitely many different values on the 
interval < t < cx). Moreover, the moment integrals 

/■oo 

lin = / t"d$(t), nGNo, (40) 
Jo 

must be positive and finite for all finite values of n. If we insert the geometric series 1/(1 -|- rr) = 
'Yl^=o{~^Y into the integral representation ([5U|) and integrate term-wise using (jlU]) . we see that a 
Stieltjes function F{z) can be expressed by its corresponding Stieltjes series: 

oo 

F{z) = ^{-IYp^z^. (41) 

u=0 
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Whether this series converges or diverges depends on the behavior of the Stieltjes moments i-in as 
n — > oo. 



If we insert the relationship 



1^=0 



1 + X 



1 + X 



(42) 



for the partial sum of the geometric series into the integral representation (|39p . we see that a 
Stieltjes function can be expressed as a partial sum of the Stieltjes series (jlT]l plus a truncation 
error term which is given as an integral (see for example Theorem 13-1 of [31]): 



Fiz) = ^(-l)V,z- + (-z)-+i / 



d$(t) 



+ zt 



argfz < vr . 



(43) 



Moreover, the truncation error term in (j43p satisfies, depending upon the value oi = arg(z), the 
following inequalities (see for example Theorem 13-2 of [31]): 



(-^) 



n+l 



CO j.n+1 



l + Zt 



< 



n+l\ 



yn+1 



sm{e) 



\e\ < tt/2, 

7r/2 <\9\<7r. 



(44) 



Notice that the expression fin+i l^""*"^! emerges as a natural remainder estimate for the partial 
sum X]I^=o i~^y °^ t'^^ input series. This observation will later be used in Sec. 12.2.31 for the 
construction of a special class of sequence transformations. Detailed discussions of the properties 
of Stieltjes series and their special role in the theory of summability can be found in Chap. 8.6 
of [20] or in Chap. 5 of [27]. 



2 Sequence Transformations 

In the current section, sequence transformations are introduced and discussed in several steps. The 
goal is both to introduce basic notions connected with the construction of convergence accelerators, 
and to enable the reader to read and understand further, more specialized literature sources. 

• In Sec. 12.11 an overview of some of the most common convergence acceleration methods is 
given, and the reader is presented with a crude quick-start guidance for the choice of suitable 
convergence acceleration methods. 

• In Sec. 12.21 selected convergence acceleration methods are discussed. Some of these profit from 
the availability of explicit remainder estimates and can sometimes lead to quite spectacular 
numerical results. 

• In Sec. 12.31 the methods are applied to a number of concrete computational problems, for 
illustrative purposes. The examples concern both straightforward applications of some of the 
methods as well as more sophisticated cases, where methods have to combined in an adaptive 
way. 

The Chapter thus provides both abstract as well as very concrete information regarding convergence 
accelerators. 
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2.1 Quick Start 



As a general preamble, it should be noticed that, as everything comes at a price, there is unfortu- 
nately no magic method working in all situations: a famous proof in Ref. [79] shows that a universal 
transformation accelerating all numerical sequences cannot exist, even if we restrict our attention 
to convergent sequences, neglecting the divergent case. 

No matter which convergence acceleration method is finally selected for a particular problem at 
hand, some practical effort is always necessary to get the best results: on the one hand, a straight- 
forward usage of the algorithms presented in this review, in the sense of computational recipes, can 
often lead to misleading and nonsensical numerical results; on the other hand, quite spectacular 
results can often be achieved if careful use is made of all the available information. A two-step 
procedure should be followed: before starting the actual numerical investigations the known theo- 
retical results should always be reviewed, to understand if by chance the problem at hand satisfies 
the applicability conditions of one or more of the existing methods. Only if the answer to this 
question is negative the user should try out some algorithms blindly, and even in this case a pre- 
liminary critical assessment of the nature of the problem can often help in immediately ruling out 
a number of methods. To illustrate these points, we lay out a possible rough guidance for choosing 
and implementing the best sequence acceleration method among those which are presented in in 
the rest of the chapter. A much more complete (and interactive) version of this decision grid can 
be found in the form of computer program named HELP as introduced and discussed in Ref. [25]. 
Here, we attempt to present a rather compact discussion. 

Before we turn our attention to specific cases, let us mention that it is possible to divide many of the 
algorithms discussed in the literature into two classes. In the case of "more traditional methods" , an 
assumption about the mathematical structure of the remainder term Vn is implicitly coded into the 
algorithm. The algorithm is constructed to eliminate remainder terms of a specific structure step 
by step. The terminology is that the sequence transformation is exact for a specific subset of input 
series. In the case of "modern nonlinear sequence transformations," an explicit remainder is coded 
directly and explicitly into the algorithm: the calculation cannot proceed without the external 
input of remainder estimates in addition to the terms of the series to be summed. In the case of 
the modern nonlinear sequence transformations, the algorithm only has to make an assumption 
regarding the correction term which relates the exact remainder and the remainder estimate, but 
it is unnecessary to eliminate the entire remainder from the series (only the correction term needs 
to be dealt with). 

For both the "modern" as well as the "traditional" algorithms discussed here, one may indicate a 
certain unifying concept. Indeed, many of the convergence acceleration methods can be cast into 
the general structure of the E algorithm (see [80,81] and Sec. 3.3 of Ref. [31]). By examination 
of the E algorithm, the following fundamental paradigm, on which all methods seem to be based, 
can be exposed in a particularly clear fashion: "Make an assumption about the remainder and 
eliminate it from the partial sums of the series as far as possible." 

In view of the "non-universality theorem" of Ref. [79], the task of finding a suitable acceleration 
transformation acquires a certain heuristic component. Instead of giving an immediate answer, we 
should therefore attempt to formulate some useful questions which should be asked first: 

• Is a remainder estimate a;„ for the remainder term r„ known, or is no such information 
available? In the first case, one has the option of using one of the nonlinear sequence trans- 
formations discussed here in Sees. 12.2.31 and 12.2.41 and in Sees. 7 and 8 of Ref. [31]. These 
transformations profit from remainder estimates and make implicit assumptions about the 
correction term which corrects the remainder estimate in the direction of the explicit re- 
mainder. If no remainder estimate is available, then one only has the option of one of the 
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more traditional methods covered in Sees. 12.2.11 and 12.2.21 which do not need as input an 
exphcit remainder estimate. Note, however, that the more traditional algorithms nevertheless 
make assumptions about the mathematical structure of the remainder term: the quality of 
the results varies strongly depending on the validity of these assumptions, and thus one can 
generally recommend to try more than one method on a particular given problem. 

• Is the convergence of the input series known to be linear or logarithmic, and is the sign pattern 
of the terms of the input series alternating or nonalternating? If this knowledge is available, 
then the choice of applicable convergence accelerators can be narrowed down further. Details 
can be found throughout Chap. [2j In general, the application of convergence accelerators to 
alternating series is numerically much more stable than to nonalternating series. Furthermore, 
the acceleration of logarithmic convergence is generally recognized as a more demanding task 
as compared to the acceleration of linear convergence. However, as discussed in this review, 
a few rather powerful algorithms are currently available. 

• Is it the aim to obtain a specific numerical result for a physical quantity under investigation, 
e.g., an energy level or a cross section, to the highest accuracy based on a finite number 
of terms of an input series, or is the goal to develop a high-performance algorithm, e.g., 
for a special function which is given by an asymptotic series expansion whose terms are 
known? In the first case, there is no need for a computationally efficient implementation 
of the acceleration algorithm: typically, a knowledge of its definition is sufficient, and in 
order to avoid numerical loss in higher transformation orders, one may use multi-precision 
arithmetic, which is implemented in today's computer algebra systems (e.g., Ref. [82]) in 
intermediate steps. E.g., if there is no concern about computational efficiency, then it is 
possible to directly key the definitions given in Eqs. (j84p and (jSSp into a computer program, 
without making use of the recursive schemes defined in Eqs. (ISip and (|82p below. The same 
applies to prototyping applications for numerical algorithms. In the second case, however, 
one typically has highly powerful recursive schemes available for the calculation of higher- 
order transformations which reduce the calculational demand by orders of magnitude. The 
latter aspect is crucial to the applicability of convergence accelerators. In order to appreciate 
this aspect, one should remember that it was the discovery of the e algorithm in Ref. [14] 
which sparked a huge interest in Fade approximants which in turn constitute a paradigmatic 
example of a versatile convergence acceleration algorithm. 

We can now turn our attention to the actual algorithms which have a certain prominent status 
among the many methods available from the literature. Here, it is inevitable to name some algo- 
rithms which will be defined only in Chap. [2] below. The authors hope that this aspect will not 
lead to confusion. In general, one may distinguish four cases, according to Eq. (f27l) : 

• Linear convergence, alternating series: In this case, a number of methods is available, e.g., the 
e algorithm or the (iterated) process always lead to some acceleration. These two algo- 
rithms are discussed in Sees. I2.2TT] and [2.2.21 respectively. A further possibility is provided by 
polynomial extrapolation with geometrically convergent interpolation points (see Sec. 12.2.21 
Sec. 6.1 of Ref. [31], and p. 6 ff. as weh as pp. 36-42 of Ref. [19], and Sec. 3.1 of Ref. [83]). 
The mentioned algorithms fall into the "traditional" category discussed above, and are by 
construction exact for model sequences which approximate typical cases of linearly conver- 
gent series. Moreover, the acceleration of alternating series typically constitutes a numerically 
very stable process in higher transformation orders. If a remainder estimate is available, then 
one may also study the behaviour of special sequence transformations, typically the nonlinear 
t and d, and the r and 6 transformations discussed here in Sees. 12.2.31 and 12.2.41 (see also 
Table [Din the current review) and in Sees. 7 and 8 of Ref. [31]. Observe, however, that the 
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more traditional algorithms also make special assumptions about the mathematical structure 
of the remainder terms, and that these may lead to numerically very favourable results if the 
series matches the corresponding assumptions of the algorithm. 

• Linear convergence, nonalternating series: From a theoretical point of view, there is no dif- 
ference between the cases of linearly convergent, alternating series and linearly convergent, 
nonalternating series. However, from a practical point of view, the difference is huge: be- 
cause any convergence accelerator automatically entails the calculation of higher weighted 
differences of elements of the input series, the calculation of higher-order transformations of 
a nonalternating series can be a numerically highly unstable process. The following question 
does not even have to be asked, as it usually penetrates practical considerations immediately: 
Does the accelerated series converge fast enough to the required accuracy in order to be able 
to overcome the pileup of numerical loss of precision due to roundoff? If the convergence 
is fast enough, then the direct application of one of the above mentioned algorithms to the 
input series is possible. If the roundoff errors pile up too fast, then it becomes necessary 
to make additional efforts. One possibility is to take the nonalternating series as input to a 
condensation transformation, to generate an alternating series as a result of the condensa- 
tion transformation, and then, to take the resulting series as input to an additional second 
acceleration process: here, variants of the CNC transformation (see Sec. 12.2.6]) may be useful. 

• Alternating marginally convergent series: If remainder estimates are available, then it is 
possible to use nonlinear sequence transformations, e.g., the u, d and v, and the y, 6 and (p 
transformations (see Sec. 12.2.41 here and Sees. 7 and 8 of Ref. [31]). Only for the alternating 
case, the e algorithm (see Sec. 12.2. ip may still lead to some acceleration, but its numerical 
results in the case of logarithmically convergent series are in general very poor [see also 
Eq. (I58D below]. 

• Logarithmic convergence, nonalternating series: This is in general the most problematic case, 
because numerical roundoff almost inevitably piles up in higher transformation orders. If it is 
possible to clearly discern the mathematical structure of the input series, then it may still be 
possible to use a one-step convergence acceleration algorithm with good effect. An example 
would be polynomial extrapolation with interpolation points which are given as inverse powers 
of the index, if the input series has a corresponding structure (see Sec. \2.2.2\ Sec. 6.1 of 
Ref. [31] and Ref. [19]). Another option is to transform first to an alternating series via a 
condensation transformation, and to use an efficient accelerator for the resulting alternating 
series. The latter method leads to variants of the CNC transformation (see Sec. I2.2.6p . 

The above list of algorithms is incomplete. E.g., the {) algorithm derived in Ref. [16,19] is known to 
be quite an efficient accelerator for both linearly convergent as well as logarithmically convergent 
sequences. In some sense, the algorithm has a special status because it cannot be derived as 
a variant of the E algorithm. As remarked in Sec. 15.5 of Ref. [31], the "& algorithm as well as 

(n) 

its iterations are less robust than the e and p algorithms, and in addition, the 'd algorithm 
is outperformed by nonlinear sequence transformations (notably, the u and v transformations) in 
many cases of logarithmically convergent sequences. Nevertheless, because of its unique position as 
a rather robust, and rather universal accelerator, the t? algorithm definitely has to be mentioned 
in this review (see Sec. I2.2.2p . In that sense, and at the risk of some over-simplification, we can 
say that the 'd algorithm interpolates between the rather universally applicable, but moderately 
powerful e and p algorithms on the one hand and the specialized, but potentially spectacularly 
powerful nonlinear sequence transformation with remainder estimates on the other hand. 

A rather heroic attempt at a universally applicable accelerator is the COMPOS algorithm presented 
in [25] , which consists of many base sequence transformations (up to 9 in the original implementa- 
tion) being performed in parallel. The results of the parallel transformations are sorted according 
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to their apparent convergence via the SELECT algorithm (here, the apparent convergence is given 
by the relative change in the numerical values of the transforms as the order of the transformation 
is increased). The output sequences (or a subset of them) are COMPOSed by making use of the 
-E-algorithm. It has been proven that if one of the base transformation accelerates a sequence, 
then also the composite transformation will accelerate the same sequence. In addition, the speed 
of convergence of the composite transformation is usually not much worse than that of the best 
base transformation for the sequence (and can even be better if many of the base transformations 
converge quickly). A not so serious criticism of the COMPOS algorithm is based on the observa- 
tion that the algorithm fails when one or more of the base transformations converge quickly to a 
wrong answer. The interested reader can consult [25] for further details, but this drawback is of 
course shared by un-composed transformations as well. A rather serious drawback of the COMPOS 
algorithm, however, lies in the computational effort implied by the necessity to perform the initial 
transformations in parallel: if the goal is to develop an efficient numerical algorithm for a math- 
ematical entity defined by a slowly convergent series by implementing a convergence accelerator, 
then the computation overhead implied by the COMPOS algorithm can be quite severe. 

Hints concerning the actual implementation of the methods discussed here can be found, e.g., in 
Refs. [31,84] concerning sequence transformations, including the e and the p algorithms, and in 
Ref. [25] concerning a number of rather sophisticated algorithms, including COMPOS. Finally, let 
us remark that a complete implementation of the CNC transform, adapted to the computation of 
the Lerch transcendent, is available for download [85]. With these preliminaries being discussed, 
let us know go in medias res. 



2.2 Selected Methods 
2.2.1 Pade Approximants 

As already discussed in Section 11.2.31 in many problems of theoretical physics, a physical quantity 
f{z) can be presented in terms of its asymptotic or perturbation expansion (|22p with respect to 
the asymptotic power series {^"'l^o- Often, the convergence problems with such power series can 
be overcome by making use of the Pade approximants. The literature on Pade approximants is 
abundant. A classic reference is [17], and further volumes (Refs. [20,27,86,87]) provide additional 
valuable information. From among the many other articles that deal with the subject, we only 
mention Refs. [20,27,88-93]. A Pade approximant produces a rational approximation to the original 
input series which is given by the ratio of two polynomials. For instance, the [l/m] Pade approximant 
to a power series (j22]) 

oo 

/(z) = 5^afez^ (45) 

fc=0 

is given by the ratio of two polynomials Pi{z) and Qm.{z) of degree I and m, respectively: 



Pijz) _ P0+P1Z + ...+ piz'- 
Qm{z) qo + qiz + ... + : 



[l/m]f{z) = = „ , „ ^ , , „ , (46) 



where qq is canonically fixed to a unit value. The polynomials Pi{z) and Qm{z) are constructed so 
that the Taylor expansion of the Pade approximant (06]) agrees with the original input series up 
to the order of Z + m 

00 

f{z) - [l/m]f{z) = J2»k^'-[l/m]f{z) = 0(z'+™+^). (47) 

fe=0 

This equation, therefore, can be used in order to find the coefficients of the polynomials Pi{z) and 
Qm{z)- That is, by comparing coefficients at the same orders of z and assuming go = 1; one gets 
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a system of / + m + 1 coupled equations to determine po,pi,qi,p2,q2, ■ ■ ■ ■ In most cases it is quite 
inconvenient to solve this system of equations. A more attractive recursive prescription for the 
evaluation of Pade approximants is the e-algorithm given by the formulas [14,31] 

eL^; = , ei") = , (48a) 

4"i = + v[er^^ - 4"^] (48b) 

(n) 

for n,k € Nq. For the case when the input data eg = s„ are the partial sums 



n 

k=0 k=0 



J2 «^ (49) 



of the power series Eq. (I22p . the elements e'^ of the two-dimensional table of transformations 
produce Pade approximants according to 



e^^J = [n + k/k]f{z). (50) 

The odd elements of the table of transformations £2^+1 (with j G Nq) are only auxiliary quantities. 

The above formulas would not be meaningful without some interpretation. First, we observe that 
the Pade approximant [n/0]f{z) exactly reproduces the partial sum Sn of the input series. For 
the computation of [n + 1/1] j(z), we need s„, Sn+i and Sn+2 (two elements in addition to the 
start element s„), and therefore the approximant [n + l/l]j(z) constitutes a transformation of 
second order. In order to compute the Pade approximant [n + k/k]f(z), we need the elements 
s„+i, . . . , Sn+2k in addition to Sn, and the transformation therefore is of order 2k. 

In Eq. (4.3-2) of Ref. [31], it is remarked that the elements of the e table connected by the formula 
(|48b|) form a rhombus in the (n, k) plane. The general paradigm of the e algorithm is to use 
the recursive scheme in order to increase the transformation order as far as possible for a given 
number of elements of the input sequence. If the elements s„, . . . , Sn+m are known, then highest- 

(n) 

order transform which can be calculated is e2[m/2|' leather sophisticated techniques are available 
for a computationally efficient implementation of recursive schemes like the one given in Eq. (j48p . 
Typically, in most cases the surprising conclusion is that a single, one-dimensional array of elements 
is sufficient in order to perform the transformation, whereas the structure of the recursion relation 
would otherwise suggest that at least a two-dimensional array is needed. Details on the latter point 
can be found, e.g., in Sees. 4.3, 7.5 and 10.4 of Ref. [31]. 

An important point, which has certainly not escaped the reader's attention, is that for z = 1, the 
two forms of the input data in Eq. ()49p become equivalent with = afe, and in this case, the Pade 
approximant is nothing but a (nonlinear) sequence transformation of the input series {s„}^q. In 
fact, the e algorithm makes no difference between input data which are partial sums of a power 
series of term akz'' and partial sums of a simple input series of terms a^. It might also be helpful 
to remark that the e algorithm is exact for sequences of the form (see e.g. Ref. [94]) 

fc-i 

Sn = ^CjX], |Ao| > |Ai| > ... > |Afc_i| , (51) 

j=0 

which is a terminating series composed of, loosely speaking, "linearly convergent terms" (if |Ao| < 
!)■ 



18 



It is time now for a very brief historical detour. The modern history of convergence started with 
the seminal paper [73], in which the Shanks transformation efc(s„) was introduced. This trans- 
formation later turned out to be equivalent to the calculation of Pade approximants. In fact, the 
transformation ek{sn) and the e algorithm are related by [see, e.g., Eq. (4.2-2) of Ref. [31]]: 

6£)=efc(s„,). (52) 

It was the article [14] which elucidated that a highly efficient recursive scheme for the calculation 
of the approximants ek{sn) can be devised. The Shanks transformation, from our point of view, 
constitutes a way of expressing a diagonal Pade approximant \n + k/k]f{z) to a function. 

The entries of the e table with odd lower index satisfy 



Let us interpret this relation. We have assumed that the sequence composed of the s„ can be 
written as {Y^^=QOikZ^ 

}"=o- The sequence composed of the As„ is thus {a„+i2;'^+^ 

)n=o- The 

expression for the nth sequence element a^+i z"^^^ can be written as the partial sum of a power 
series, 



n 



Asn = ao + Y^ (afc+i - ak z^j . (54) 

A;=0 

As n — > oo, all coefficients of this power series vanish, and the function Asqo is zero. The quantities 
ej(Asn) can thus be identified as [n -|- j'/j]-Pade approximants to a specific function defined by a 
power series whose partial sums are given by Eq. (j54p . and since all coefficients of that function 
vanish as n — > oo, we have ej(As„) = 0{z'^~^'^^^^) by construction of the Pade approximant. 

In writing expressions like e"2k ^k{sn), we encounter a certain problem in the notation. Indeed, 
n is an index which is normally associated to a specific element s„,, and once n is fixed, s„, is fixed 
also. Here, as in many other literature sources in the field, n is often used in a double meaning: 
first, to indicate the starting element of a transformation (fixed index), and second, to indicate 
all the sequence elements Sm with m G {n, . . . , n + j} which are needed in order to construct 
a transformation of a specified order j. Observe that in order to calculate ek{sn), the elements 
Sn, ■ ■ ■ , Sn+2k are needed, as ek{sn) is a transformation of order 2k. The simple notation ek{sn) might 
otherwise suggest that the transformation takes as input only the element s„, and this is simply 
not the case. However, a certain compromise between the exactness and completeness of a notation 
and its compactness and practical usefulness must be found, and the notation ek{sn, ■ ■ ■ , Sn+2k) 
undeniably looks rather clumsy. The customary notation efc(s„,) is thus being adopted here. 

Let us briefly review the behaviour of the e algorithm for a number of model series. We follow 
Ref. [94]. For a linearly convergent series with remainders 



+E 



]cj Xj , 1 > Ao > Ai > • • • > , n ^ oo , (55) 
the e algorithm has the following behaviour as n ^ oo for fixed k, 

^2k ~ s + Cfc -ifZTT: ^fc ' n ^ oo . (56) 



The e algorithm thus accelerates convergence according to the definition in Eq. (j35|) because 
(Afc/Ao)"" — > for n — > oo (we have A^ < Ao by assumption). For a marginally convergent se- 
quence with strictly alternating remainders. 



oo 



« + (-i)"Z.(^Tt)^' ^^^^ 
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the following estimate has been found in Ref. [94] for fixed k, assuming cq 7^ 0, 



(„) (-l)"(fc!)^ .... 

4fc + 22fc (/3 + ^)2fc+i ^o- ^^oo- (58) 

Also here, the e algorithm accelerates convergence because (/3 + n)/(/? + n)^'^"'"^ ^ for n ^ 00 
for k > 0. However, the expression (A;!)^ in the numerator suggests that the e algorithm encounters 
problems when one tries to increase the order of the transformation k for given index n of the start 
element of the transformation. For monotone (nonalternating) series, 

00 

Sn^ s + ^ I ' /3>0, n^oo, (59) 

^ + n)J+i 

the e algorithm is unable to accelerate convergence. 



2'^ (A; + 1) (/3 + n) ' ^ ^ 



Indeed, the convergence theory of Fade approximants is quite well developed [17,27], and the above 
examples are provided only for illustrative purposes. 

Besides the model series discussed above, we would like to recall briefly the convergence and accel- 
eration results obtained with the e algorithm when applied to two important classes of sequences: 
totally monotonic and totally oscillating sequences. Here, we start our discussion with totally 
monotonic sequences Sn which by definition have the following property: 

(_l)fcA^s„>0 VA;,nGNo. (61) 

This condition implies that since s„ > (set /c = in Eq. (I6ip ). we have Sn > Sn+i > (set k = 1 
in Eq. (|6ip ). A totally monotonic sequence converges in such a way that Vn, we have s„ > s > 0. 
Another, equivalent definition of a totally monotonic sequence is possible. Namely, it can be proven 
that a necessary and sufficient condition that s„ is a totally monotonic sequence is given by the 
existence of a bounded and non-decreasing measure da(x) in [0, 1] such that: 



x"da(a;) . (62) 







The sequences Sn — A*^ where < A < 1 and — l/(^ ~l~ 1) '^^^ examples of totally monotonic 
sequences. 

The properties of totally monotonic sequences have been discussed in detail elsewhere (see, for 
example, Refs. [25,95]). In the present manuscript, we just mention the convergence properties. 
Namely, it has been proven in [96] that if the e algorithm is applied to a sequence s„ which converges 
to s and if there exist parameters a and b such that a s„ + 6 is a totally monotonic sequence, then it 
is possible to prove the following properties which ensure the convergence of the Fade approximants 
given by the e algorithm to the limit s, as both n and/or k are increased. Namely, we have Vfc and 
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Vn (see p. 89 of Ref. [25]), 



< 


«e2fc+2 + ^ < «e2fc + b , 


(63a) 


< 




(63b) 


< 


«4fc+2 + ^^"4fc 


(63c) 


< 


(n) , I ^ (n+2) , , 


(63d) 


VA; 


fixed, lim el?^ = s , 

n— >co '"^ 


(63e) 


Vn 


fixed, lim = s . 

fe— »oo 


(63f) 



As expressed by Theorem 2.23 of Ref. [25], all the columns and all the descending diagonals of the 
e algorithm can be used for the acceleration of convergence of a totally monotonic series, provided 
the input series is a non-logarithmic totally monotonic series. We recall that the e algorithm is 
unable to accelerate logarithmic convergence. 

In a similar way as for a totally monotonic series, the e algorithm can also be useful for the 
summation of totally oscillating series. Sequences of that kind are defined in the following way: 
a Sn is totally oscillating if the sequence (— l)"'s„ is totally monotonic. Again, it has been proven 
that if the e algorithm is applied to a sequence Sn which converges to s and if there exists a and b 
such that asn + 6 is totally oscillating sequence, then we have V/c and Vn according to Ref. [25]: 



0<«4f+2 + ^<«4r + ^> (64a) 

a e£"+') + b<a ef^,'^ + 6 < , (64b) 

«(4r^^-4ro<«(4i;v^-4S2)<o, (64c) 

< a (4l;?^ - < a {4r' - e?r^) , (64d) 

0<aeg% + b<aeg^+'^ + b, (64e) 

a e£"+') + b<a eg"/^'^ + 6 < , (64f ) 

(n) 

\/k fixed, lim = s , (64g) 

Vn fixed, lim e^2k ~ ^ ■ (64h) 



2.2.2 Further Convergence Accelerators w^ithout Remainder Estimates 

The e algorithm takes as input only the elements of the series to be accelerated. No additional 
remainder estimates are needed. This property is shared by the process, the p algorithm, and 
the "i? algorithm, which are discussed in the current subsection. 
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The process is widely applied to the summation of slowly convergent series. It is usually at- 
tributed to Aitken [4], but there are indications that it may already have been known to Kum- 
mer [97]. Within this process, one makes use of the forward difference operator (jlip in order to 
transform the original sequence {sn}^=Q into a new sequence {An}^=o, whose elements are defined 
by: 

An = Sn- . (65) 

The structure of this equation clearly explains why Aitken's method is called method. Before 
discussing the convergence properties of the transformation (|65p. we would like to mention that 
originally Aitken's process was designed to deal with the model sequence of the following form [4], 
which is a truncated form of Eq. ()5ip : 



s„ = s + cA". (66) 

As seen from this equation, each sequence element s„ contains two coefficients c 7^ and |A| < 1 
as well as the limit s. The first and the second forward differences of such a sequence read: 

As„ = cA" (A - 1) , A^s^ = cA" (A - 1)^ . (67) 

By inserting now these relations into Eq. (j65p . we may easily find: 

An = (s + cA) — cA = s, 



which implies that the A^ process is exact for the model sequence ()66p at every stage. Like any 
other algorithm of numerical analysis, this method has its own domain of validity. Because the 
(convergence) properties of the transformation ()65p have been discussed in detail elsewhere [10, 18, 
19], here we just mention that the A^ process (i) accelerates linear convergence and (ii) is regular 
but not accelerative for logarithmically convergent sequences. This indicates that A^ process has 
similar properties as the e-algorithm but is usually not as powerful as the latter one [31]. One 
may increase the accelerative power of the Aitken's method, however, by applying this method 
iteratively. That is, it is first applied to the original sequence {sn}^0' then to the transformed 
sequence {^„}^q and so on. The detailed discussion on the iterated A^ process is given in the 
texts by Brezinski and Redivo Zaglia (see pp. 128-131 of Ref. [25]) and Weniger (see pp. 223- 
228 of Ref. [31]) and, therefore, will be omitted in the present report. Apart from the iterative 
procedure, moreover, a large number of modifications to the A^ algorithm have been presented 
in Refs. [98-100] and many others. At the end of a rather brief discussion on the Aitken's A^ 
method, we shall note that this method has very close connection to the e algorithm as well as to 
the Fade approximants. For instance, it immediately follows from the derivation of the sequence 
transformation An via the model sequence ([66|) that An = 62"^ = [n + l/l]j(2;) is effectively a Fade 
approximant to the input sequence given by s„ = ^11=0 (^kZ^ ■ 

Beside the A^ process, another very well-known sequence transformation is the p algorithm [15], 
which is exact for logarithmically convergent sequences of the following structure. 



sJl + aix^ ^ + ... + ak 
xi + hixi'^ + ... + bk 



/c,neNo, (69) 



where the x„ are usually chosen to be 2;„ = /3 + n with /? > 0. The p algorithm is defined by the 
scheme 

pLI = 0, pS") = sn^ (70a) 

(n) _ (n+1) , Xn+k+1 - Xn , . 

Pk+i - Pk^i + („+i) („) ■ y'^") 

Pk -Pk 
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The recursive schemes for the e and p methods are almost identical, as seen from Eqs. (j48p and (j70p . 
The only difference is that p algorithm also involves a sequence of interpolation points {x„}^q, 
which are usually chosen to satisfy the condition < xq < xi < 3:2.... and x„ — > 00 if n — > 00 (see 
p. 231 of Ref. [31]). Moreover, similar to the case of the e algorithm, only the elements with 

(n) 

even orders serve as approximations to the limit while the elements P2fc+i ^'^^^ '^^^ orders are just 
auxiliary quantities which diverge if the whole process converges. Despite their formal similarity, 
the e and p methods differ significantly in their ability of convergence acceleration. For example, 
the e algorithm is not able to accelerate nonalternating logarithmically convergent sequences [see 
Eqs. ([59]) and <^ and the discussion around Eq. (4.2-8) and (4.2-9) of Ref. [31], as well as Ref. [94]]. 
In contrast, such a logarithmic convergence can be accelerated by the p algorithm, which fails in 
the case of linear convergence. Thus, the two algorithms can be considered to some extent as 
complement ary. 

A scheme which interpolates between the e and p algorithms is the "Q algorithm [80,81]. It assumes a 
special place in the hierarchy of the algorithms as it cannot be derived on the basis of the otherwise 
very general E algorithm discussed below in Sec. 12.2.71 The ??-algorithm is defined by the scheme 

^g=0, = (71a) 
= 4^1^ + -^-^ TT > (71b) 

2k+l 2k-l „{n+l) _ „(n) ' ^ ^ 



(„oin+2) _ nin+l)\ ( a(n+2) _ .(n+l)\ 
An) _ An+1) ^2k ) y^2k+l ^2k+l ) 

^2k+2-^2k + fAn+2) _^An+l) ..(n) ^ " ^^'^^ 

y"2k+l "2k^\ "T '^2k+l^ 

with k,n £ Nq. Here, we have written out all expressions which might otherwise be formulated a 
little more compactly by using the difference operator A. The algorithm can be considered as 
an intermediate form of the e and the p algorithm, as it is able to accelerate both some linearly 
and some logarithmically convergent sequences. Although the e-algorithm and the p-algorithm are 
particular cases of the £^-algorithm, the i?-algorithm is not, and the same can be said about various 
sequence transformations which were obtained by a manual modification or iterations of some other 
previously known transformations, e.g., of the 1) algorithm. For a complete list of such sequences 
we refer the interested reader to [25,31]. 



2.2.3 Sequence Transformations with Remainder Estimates 

As seen from Eq. ([16|) . the coefficients ao,ai, . . . ,ai+m of the formal power series (p2]) suffice 
to determine the Fade approximant [l/m]f{z) completely without the need for any additional 
information. The same can be said about the methods discussed in Sec. 12.2.21 all the input required 
is in the partial sums of the input series. However, the estimate of a truncation error (remainder 
term) can often be given as the first term of asymptotic power series not included in the partial 
sum (see Sec. 11.2.6) and Refs. [76,101]). In this Section, we discuss the sequence transformations 
which can benefit from such an additional information via remainder estimates. One of these 
transformations gave very good results in the case of the sextic anharmonic oscillator, and there 
is strong numerical evidence that it is able to sum the extremely violently divergent perturbation 
series for the octic anharmonic oscillator [102]. Nonlinear sequence transformations have been the 
subject ofanumber of monographs [1,10,18,19,21,22,24-26,37] and review articles [30,31,103,104]. 

In order to start the discussion on the sequence transformations with explicit remainder estimates, 
let us again consider the sequence {sn}^o °^ partial sums of an infinite series either converges to 
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some limit s, or, if it diverges, can be summed by an appropriate summation method to give the 
generahzed hmit s. As seen from Eq. (I23p . the summation of slowly convergent or divergent series 
to the (generalized) limit s requires the estimation of the remainders r„ from the sequence elements 
Sn- Unfortunately, a straightforward elimination of the remainders is normally not possible, since 
they are either unknown or not easily accessible. Thus, in realistic situations one can only hope 
to eliminate approximations to the actual remainders. If this can be accomplished, {s„}^q is 
transformed into a new sequence {s^}^q, whose elements can for all n > be partitioned into the 
generalized limit s of the original sequence and a transformed remainder r^: 

s'n = s + r'^. (72) 

This approximative elimination process was successful if the transformed remainders r'^ vanish as 
n oo, which implies that the transformed sequence {s^}^o converges to s. Thus, a sequence 
transformation is a rule which transforms a slowly convergent or divergent sequence {s„}J^q into 
a new sequence {s'^j^g with hopefully better numerical properties. 

The question arises how assumptions about the structure of the remainder estimate can be trans- 
lated into an efficient computational scheme. The direct elimination of approximations to r„ from 
Sn can be very difficult. A considerable simplification can often be achieved by means of a suitable 
reformulation. Let us consider the following model sequence (see Section 3.2 of [31]): 

Sn = S +rn = S + UJnZn- (73) 

Here, cOn is a remainder estimate [cf. Eq. (|26p ]. which has to be chosen according to some rule and 
which may depend on n in a very complicated way, and Zn is a correction term, which should be 
chosen in such a way that it depends on n in a relatively smooth way. Moreover, the products uJnZn 
should be capable of producing sufficiently accurate approximations to the actual remainders 
of the sequence to be transformed. 

The ansatz (|73p has another advantage: One may now easily devise a systematic way of construct- 
ing a sequence transformation which is exact for this model sequence, if one can find a linear 
operator T which annihilates the correction term z„. Just apply T to — s]/u;„ = Since 
T annihilates Zn according to T(zn) = and is by assumption linear, we find that the following 
sequence transformation T is exact for the model sequence (f73]) . according to Eq. (3.2-11) of [31]: 

T{Sn,UJn) = = S. (74) 

T{l/uJn) 

Simple and yet very powerful sequence transformations are obtained if the annihilation operators 
are based upon the finite difference operator A defined in Eq. pip . As is well known, a polynomial 
Pk-ii'iT') of degree A; — 1 in n is annihilated by the k-th power of A. Thus, we now assume that the 
correction terms {znj^o chosen in such a way that multiplication of Zn by some suitable 

quantity Wk{n) yields a polynomial Pfc_i(n) of degree /c — 1 in n: 

Wk{n) Zn = -Pfc-i(ra) . (75) 

Since A'^Pfc_i(n) = 0, the weighted difference operator T = A^Wk{n) annihilates Zn, and the 
corresponding sequence transformation (j74p is given by the ratio 

n-Hf ( M \ A^{wfc(n)s„/u;„} 

T; (...(n)|.„,.„) = ^,^^^^^^1^^^ ■ (76) 

Indeed, Eq. (|76p represents the general form of the nonlinear sequence transformation, based on 
annihilation operator T. The exact form of this transformation depends on the particular choice of 
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the Wk{n). For instance, Wk{n) = {n + (3)^ ^, which annihilates correction terms given by inverse 
power series, 

with /3 > 0, yields the following sequence transformation [31,105]: 



^ " AM(n + /9)^-V^n} 

^ (_l)i ^ ^ (^ + n + j)^"^ s„+j 



j=0 



- , (78) 



^ ^ k \ (/3 + n + i)'=-i 



where the subscript k denotes the order of the transformation, and n denotes the starting index 
of the input data (fc + 1 sequence elements s^,, . . . Sn-\-k 

and remainder estimates a;„ . . . uJn+k are 
needed for the computation of the transformation). The shift parameter /3 has to be positive in 
order to admit n = in ()78p . The most obvious choice would be f3 = 1. 

Apart from the transformation (j78p . we may also adopt the unspecified weights Wk{n) in (j76p to 
be the Pochhammer symbols according to Wk{n) = (n + f3)k-i = r(n + (3 + k — l)/r(n + /?) with 
/3 > 0. This transformation is exact for corrections terms of the following structure, 

-|(^- <™> 

This yields the transformation (see Eq. (8.2-7) of [31]): 

Mn),a ^ _ A^' {(n + ^)fc_i gn/^n} 

A'^{(n + /?)fc-l/t^n} 



k 
j=0 



k \ i(3 + n + j)k 



k 
3=0 



j J (/3 + n + /i;)fc_i Wn+i 

(80) 



k \ {P + n + j)k-i 1 



j J {P + n + k)k-i uJn+j 



Again, the most obvious choice for the shift parameter is /3 = 1. The numerator and denominator 
sums in (j78p and (jSOp can be computed more effectively with the help of the three-term recursions, 
according to Eq. (7.2-8) of [31], 

Vi(/3) - ifc (/?) (;3 + n + A, + l)fc 

and according to Eq. (8.3-7) of [31], 

cW - (/3 + n + fc)(/3 + n + fc- 1) c^(n).^. . ^ 

S,+i(/3) - S, (/?) - (^^^^2fc)(/3 + n + 2fc-l) ' ^^^^ 

The initial values Lq"''(/3) = S'q"^(/3) = s„/a;„ and Lq"'*(/3) = sl^\f3) = l/a;„ produce the numerator 
and denominator sums, respectively, of C^i^\p, Sn,iOn) and sj^\p, Sn,u)n)- 
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The remainder estimate 

LOn = an+l = ^Sn, (83) 

was proposed by Smith and Ford in Ref. [106] for both convergent as well as divergent series with 
strictly alternating terms. The use of this remainder estimate in ()78p and ()80p yields the following 
sequence transformations, defined according to Eqs. (7.3-9) and (8.4-4) of [31]: 



J2{-iy ( ^ ^ {fi + n + jf-^ 
i=o 



3 J (/? + n + A;)'=-i Asn+j 



k 

E(-i) 

3=0 



k\ (/? + n + j) 



k-l 



(84) 



j J ip + n + k)''-^ Asn+j 



and 



k \ {l3 + n + j)k — 1 Sn+j 



k 



j=0 



j J {l3 + n + k)k-i Asn+j 



k 



j=0 



k \ {(5 + n + j)k-i 1 



(85) 



3 J {l3 + n + k)k^i Asn+j 



In Appendix A of [34], it was explained why the transformations d^"^ and (5^"^ are more effective 
in the case of alternating series than in the case of nonalternating series. The failure of the trans- 
formation for cases where the remainder estimates are invalid is also illustrated in Table 8.7 of 
Ref. [107]. Other sequence transformations, which are also special cases of the general sequence 
transformation ([76]), can be found in Sections 7 — 9 of [31] and in Section 2.7 of Ref. [25]. A sequence 
transformation, which interpolates between the sequence transformations (|78p and ()80p . has been 
derived in [108]. 



2.2.4 Further Sequence Transformations 

According to Eq. (7.3-13) of Ref. [31] (see also p. 19 of Ref. [10] and the elucidating discussion in 
Ref. [109]), one can heuristically motivate different remainder estimate for specific types of input 
series. Let us begin with series Sn = X^^Lo where the terms have the following structure, 

o-n ^ A" n® , n — > oo . (86) 

Here, oq is a constant. For |A| < 1, which can be seen as a paradigmatic case of linear convergence, 
we have in leading order 

Tn ~ ao ^ _ , n ^ oo . (87) 

Hence, we have in this case 

A 

Tn ~ an, n ^ oo . 

The prefactor is irrelevant because it is n-independent, and therefore the choice 
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for a simple remainder estimate is heuristically motivated for linearly convergent series. This choice 
leads to the t and the r transformations [see Eqs. (7.3-7) and (8.4-3) in [31] and Table [1] in this 
review] . 

We now turn to the case 

a„ ~ Qo , (90) 

with Re0 < —1, which can be seen as a paradigmatic case of logarithmic convergence. The 
remainders fulfill 

„e+i 

rn ~ ao Q-jp[ ' Qo. (91) 

This can be related to the follows, 

1 , , 

'^n'^Q-^nan n^oo. (92) 

Both remainders have to be negative if the input series consists only of positive terms, assuming 
the series terms to be real. This follows from the definition Sn = s + rn, which can be reformulated 
as r„ = —YlT=n+i^k for s„, = J2k=o^k- For the case ([9T]) . the negativity in ensured because 
of the condition Re0 = < — 1 which is naturally imposed if the input series is assumed to 
be convergent. The prefactor is irrelevant because it is n- independent, and therefore the choice 
LOn = nan, or a slight generalization, 

uJn = {P + n) an , (93) 

is heuristically motivated for logarithmically convergent, nonalternating series. By shifting n — > 
n + (5, we have avoided potential problems in the non-asymptotic region (notably, the case n = 0). 
The choice (I93p leads to the u and the y transformations [see Table [Hand Eqs. (7.3-5) and (8.4-2) 
of Ref. [31]]. 

One may refine this choice a little bit, by observing that for an ~ a^nP , we have 

"'^'^'^^^ -^nan. (94) 



The n-independent prefactor is again completely irrelevant, and the choice 

iOn = (95) 

leads to the v and ip transformations [see Table [1] and Eqs. (7.3-11) and (8.4-5) of Ref. [31]]. 
The averaging over n implied by the more sophisticated choice (j95p as compared to (193^ can be 
advantageous. 

We now turn to the case 

a„ ~ao(-l)"n®, (96) 

with ReG < —1, which can be seen as a paradigmatic case of alternating marginal convergence 
(see also Sect. [L23]) . The remainders fulfill 

rn~ ^ao(-l)" (l + n)®, n ^ cx) . (97) 
This can be related to the follows, 

1 , , 

n^oo. (98) 
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Table 1: Variants of the transformations (j78p and (jSOp with different remainder estimates 
ujn- The assumed structure of the correction term is either an inverse power series, ac- 
cording to Eq. ([77|) . or an inverse factorial series, according to Eq. ([79]) . The main type of 
convergence corresponds to a typical model series by which the choice of the remainder es- 
timate is motivated. However, according to Ref. [31], the u, y, v and ip transformations have 
also been applied to linearly convergent series with numerically favourable results. The d 
and 6 transformations are constructed according to an alternating marginally convergent 
model series. However, the remainder estimate employed is also applicable to alternating 
linearly convergent series, and even to factorially divergent, alternating Stieltjes series, 
following the discussion in Sec. 11.2.61 [see also Eq. (j44p ]. The c and x transformations are 
proposed here as examples of alternative transformations which can easily be obtained as 
generalizations of other, existing methods. 



Assumed correction term: 
Remainder estimate Inverse powers Inverse factorials Main type of convergence 



u 



logarithmic 
(linear) 



UJr, 



0-n O-n+l 
0-n — 



logarithmic 
(linear) 



— CLn 



t 

d 



T 

5 



linear 

alternating marginal 
(linear, factorially divergent) 



OJr. 



0"n Qn+2 



X 



alternating marginal 
(linear, factorially divergent) 



The prefactor is irrelevant because it is n- independent, and therefore the choice 

tdn = a-n+i , (99) 

is heuristically motivated for logarithmically convergent, nonalternating series. The choice ()99h 
leads to the d and the 5 transformations [again, see Table [1] and Eqs. (7.3-9) as well as (8.4-4) of 
Ref. [31]]. 

The following "averaged" choice which generalizes the estimate u;^ = o-n+i) 

O'nan+2 ,^ „„N 

= , (100) 

leads to transformations which we would like to refer to as c and 



cl'^^ (/?,.„) = ^ -— — , (101) 

V f-lV 1^1 + ^ + an+j+i 
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and 



k \ iP + n + j)k 



k 

y i-iy , 

, , \ j J {P + n + k)k-i a-n+i an+j+2 

X^:\P,Sn) = ^ — ^ . (102) 



j=0 



k \ (/3 + n + j)fc„i On+j+i 



j J {P + n + k)k^i an+jan+j+2 



Not much has to be said about the c and x transformations since in our numerical studies we found 
that they appear to have very similar properties as compared to the d and 6 transformations [see 
Eqs. (j84p and (jSSp . respectively]. Other recent developments with regard to the construction of 
novel sequence transformations are described in Refs. [109,110]. As possible further modifications 
of the nonlinear sequence transformations, we should also mention multi-point methods (see, e.g., 
Ref. [111]). 



2.2.5 Accuracy Through Order Relations 

We would now like to insert a small detour into the analysis of the sequence transformation, by 
discussing so-called accuracy-through-order relations fulfilled by the rational approximants. On the 
basis of these relations, predictions for unknown series coefficients can be made. 

According to (f47|) . the difference between f{z), defined according to Eq. (f22]) . and the Fade approx- 
imant [l/m]f{z) is of order 0(z'+'""'"^) as z — > 0. This means that the Taylor expansion of [l/'m\f{z) 
around z = reproduces the first l + m + 1 coefficients oq, ai, . . . , ai+m of the perturbation series 
([221) . Consequently, the Fade approximant can be expressed as follows: 

l+m 

[l/m]f{z) = ^ Ukz'' + TZi+m+i{z) = fm+l{z) + TZi+rn+l{z), (103) 
fc=0 

where the truncation error term 7^;_|_m+i(-2) is of order 0(2:'"'"'""'"^) as z ^ 0. By performing now a 
formal power series expansion of this error term around z = 0: 

7^;+^+l (z) = 2'+"^+^ {ai+m+1 + z + . . .} . (104) 

we may find those coefficients of the perturbation series (j22p that were not used for the construction 
of the Fade approximant \l/m]f[z). 

Moreover, as it follows from the definition of the asymptotic sequence in sense of Foincare, the 
/(z) can be expressed as the partial sum fn+i{z) of the perturbation series (j22]) plus an error term 
that is also of order ©(z'^'""'"^) as z — > oo: 

f{z) = y^anZ"^ +Ri+rn+l{z) = fra+l{z) + i?i+m+l(^)- (105) 
n=0 

Similarly to the Fade approximant, the formal Taylor expansion of this error term gives: 

Ri+m.+i{z) = z'+™+^ {ai+m+1 + ai+m+2 z + . . .} . (106) 

Let us now assume that [l/m]f{z) approximates f{z) for sufficiently large values of I and m with 
satisfactory accuracy, 

f{z) ^ [l/m]f{z). (107) 
By making use of Eqs. (I103P and (I105p . this expression implies that 

Rl+rn+l{z) ~ Tll+ni+l{z), (108) 
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and, therefore: 

az+m+l ~ "Z+m+l ) "Z+m+2 ~ Otl+rn+2 , ■ ■ ■ ■ (109) 

As seen from these relations, at least the leading coefficients of the expansion for TZi+rn+i{z) 
should provide reasonably accurate approximations to the analogous coefficients of the expansion 

for Ri+,rn+l{z). 

The set of equations (I109p represents the order-through-accuracy relations for the case of Pade 
approximant. Similar relations can be obtained for all transformations listed in Table [TJ For the d 
and 5 transformations, the final results read [112]: 

f{z) - = 0(/+"+2). (110) 

The following example illustrates the order-through- accuracy relations. First, we consider an al- 
ternating asymptotic series 

oo 

A{z) ~ ^(-l)"n!z", (111) 

n=0 

whose first terms read 

A{z) = l-z + 2z^-6z^+ 24 +0{z^) . (112) 

exact term 

The second-order 6 transformation is given by 

5f = , ■ (113) 

Upon re-expansion about z = we obtain 

^ i_z^2z^-6z^+ 20/ +0{z^). (114) 

by reexpansion 

The series expansion for A{z) up to the order of z^ is reproduced (A; -|- 1 = 3), and the term of 
order z^ has to be interpreted as a prediction for the term in this order of perturbation theory. In 
higher order of transformation, the predictions become more and more accurate. 

We now consider the nonalternating series 

oo 



n=0 

which is generated by expansion of a certain integral, see Eq. (jl79p below. We attempt to illustrate 
that as long as the coupling is held variable, the prediction does not depend on the sign pattern of 
the coefficients. Namely, we have 

J\f{z) = l + z + 2z^ + 6z^ + 2Az^ + 0{z^). (116) 

The second-order S transformation is given by 



(0) _ l-3z 
^2 - i_4^ + 2z2- ^117^ 



Upon re-expansion about z = we obtain 

4°^ = i + z + 2z^ + 6z^ + 20z^ + O{z^), (118) 
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where again the coefficients of are seen to differ by 20%. The prediction of perturbative coef- 
ficients may be possible even in cases where the resummation of the series fails due to its nonal- 
ternating character (see also Table 8.7 of Ref. [107]). For nonalternating series, the prediction of 
higher-order coefficients may be possible, but the terms associated with these predictions cannot be 
used as remainder estimates. These observations have been used for the prediction of higher-order 
perturbative coefficients [113-143]. 

2.2.6 Combined Nonlinear Condensation Transformations 

Quite recently, the combined nonlinear-condensation transformation [34] (CNCT) has been pro- 
posed as a computational tool for the numerical evaluation of slowly convergent, nonalternating 
series. The idea is to divide the acceleration process in two steps. The first step, which is a re- 
ordering process, consists in a rearrangement of the terms of the input series into an alternating 
series via a condensation transformation [144] (the condensation transformation was apparently 
first mentioned on p. 126 of Ref. [145] and published only later [144]). In typical cases, the output 
of the first step is an alternating series whose convergence is marginal (i.e. p = -1). It could ap- 
pear that nothing substantial has been achieved in the first step of the CNCT. The first step of 
the CNCT merely represents a "computational investment" with the intention of transforming the 
nonalternating input series into a form which is more amenable to the acceleration of convergence. 
The second step, which represents a convergence acceleration process, consists in the application 
of a powerful nonlinear sequence transformation for the acceleration of the convergence of the 
alternating series which results from the first step of the CNCT. 

Following Ref. [144], we transform the nonalternating input series [we have a{k) = a^] 

oo 

^a{k), a(A;)>0, (119) 

fc=0 

whose partial sums are given by ([4]), into an alternating series Yl'^=Q{~^y -^j- After the first step 
of the transformation, the limit of the input series is recovered according to 

oo oo 

^o(fc) = ^(_l).A,. (120) 

fc=0 j=0 

The quantities Aj are defined according to 

oo 
fc=0 

where 

hf = 2'' a(2^ (i + 1) - 1) . (122) 

Obviously, the terms Aj defined in Eq. (jl2ip are all positive if the terms a{k) of the original series 
are all positive. The Aj are referred to as the condensed series [34], and the series Yl'j^oi~^y -^j 
is referred to as the transformed alternating series, or alternatively as the Van Wijngaarden trans- 
formed series. 

The summation over k in Eq. (I12ip dose not pose numerical problems. Specifically, it can be easily 
shown in many cases of practical importance that the convergence of ^1^=0 ^) linear 

even if the convergence of Yl^=o ^(^) only logarithmic. We will illustrate this statement by way 
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(121) 



of two examples. As a first example, we consider a logarithmically convergent input series whose 
terms behave asymptotically for A; — > oo as a{k) ~ k~^~'' with e > 0. In this case, the partial sums 

n 

a(.") = V bl^) (123) 



k=0 



converge linearly with 



a("+^)-a- 1 

lim = ^ < 1 , a{k) ~ k-'~^ , A; ^ oo. (124) 

J •' 

As a second example, we can give a series with a{k) ~ k^r'' where < r < 1 and /? real. Here, we 
have p = r < 1, and the series is (formally) linearly convergent. However, slow convergence may 
result if p is close to one. In this case, the condensed series are very rapidly convergent, 

lim = 0, a(k) ~ k^r'' , k ^ oo. (125) 

Therefore, when summing over k in evaluating the condensed series according to Eq. ()12ip . it is in 
many cases sufficient to sum the condensed series term by term, and no further acceleration of the 
convergence is required. 

As shown in [34, 146], the condensation transformation defined according to Eqs. (I120p - (ll22p is 
essentially a reordering of the terms of the input series X]fcLo^(^)' Furthermore, according to 
the Corollary on p. 92 of Ref. [146], for nonalternating convergent series whose terms decrease in 
magnitude (|a(fc)| > \a{k + 1)|), the equality (|120p holds. 

Note that the property 

A2,_i = ^ [Aj_i - a{j - 1)] , J e N , (126) 

derived in [146], facilitates the numerical evaluation of a set of condensed series, by reducing the 
evaluation of condensed series of odd index to a trivial computation. 

In the second step of the CNCT, the convergence of the series which results from the condensation 
transformation, 

oo 

^(-l)^A,, (127) 

j=0 

and which on the right-hand side of Eq. (I120p , is accelerated by a suitable nonlinear sequence trans- 
formation, e.g., by the 5 transformation (j85p . The remainder estimates are constructed according 
to 

s„ ^ S„ , ujn = As„ ^ AS„ = (-1)^+1 A„+i . (128) 

By making — if necessary — an appropriate shift of the index of the partial sums, it is possible to 
always assume that So is the initial element of the second transformation. This leads to the CNC 
transforms 

'rcNc(n) = 5f (l,So) . (129) 
These require as input the elements {So, . . . , S„, Sri+i} of the condensation transformed series. 

We here also indicate that according to Table [H it is possible to immediately generalize the CNC- 
type transformations to the the CNC-ii, CNC-y and the CNC-u, CNC-(/? transformations. The 
CNC-d and CNC-5 have been studied in Refs. [34, 147, 148]. 
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Apart from the CNCT, as discussed here, a number of other summation techniques are known which 
also make use of the division of the acceleration process into two steps. For instance, the summation 
of logarithmic sequences can sometimes be accelerated by applying an "extraction procedure." 
The main idea of this procedure is to extract from a given logarithmic sequence a subsequence 
converging linearly (on the level of the partial sums s„, not on the level of the series terms o^), 
and then to accelerate the sequence {sn}^=Q by accelerating the extracted, linearly convergent 
subsequence {s;v'(n)}5?Lo) where the M{n) are the indices of the elements of the subsequence. A 
detailed discussion of this method is given, for example in Refs. [25,149]. In particular, it has been 
proven that the acceleration of the extracted subsequence implies the acceleration of the initial 
(monotone logarithmic) sequence. Based on this result, a few algorithms have been proposed for 
the extraction of (linearly convergent) subsequences, i.e. for the determination of the N{n). Here, 
we present one of those algorithms whose aim is to find a linearly convergent subsequence where 
the parameter p is understood as defined in Eq. 

1. Let us choose 0< p < 1, and set AA(0) = 0, n = 0. 

2. Increase n — > n + 1, and set N{n + 1) = N{n) + 1, then go to step 4. 

3. Increase J\f{n) — > J\f{n) + 1, then go to step 4. 

4. If |s;v'(n)+i ~ ^N{n)\ < -^"I'Si — Sol, go to Step 2, otherwise go to step 3. 

For this algorithm, if applied to a monotone logarithmic sequence with lim„^oo Sn = s and 
lim„_+oo ^^n+i/Asn = 1, the resulting subsequence is found to satisfy [25]: 

hm ^A^(n+l)+l-^^(»+l) (130) 

Moreover, although a general proof has not been given in the general case, for many examples 
considered in [106] the following relation, which is a little stronger than (jl30p . was obtained: 

lin, f^^(!!±i)^ = A < 1 , (131) 

which immediately implies the the extracted subsequence is linearly convergent. To accelerate the 
convergence of the extracted sequence, various standard algorithms as e.g. Aitken's process can 
be applied. 



2.2.7 General E Algorithm 



A beautiful mathematical construction which deserves a particular mention is that of the E algo- 
rithm [80,81]. Rather than a traditional sequence transformation, the E algorithm is a template 
acceleration algorithm built in such a way to be exact for all sequences of the form 

Sn = s + aigi{n) + . . . + ak Qk (n) , (132) 

where the gi^s are given auxiliary sequences: with different choices for the gi^s, different acceleration 
algorithms are obtained. In particular, many of the best known sequence transformations can be 
recovered as particular cases of the E algorithm. The purpose of the E algorithm is to successively 
eliminate the gi terms in order to recover the (anti-)limit s from the 57^. 

It is interesting to examine in some detail the derivation of the E algorithm, since the employed 
method is quite general, and can be used by the reader to derive new sequence transformations 
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adapted to specific problems. A more comprehensive description can be found for example in 
Refs. [25,80]. The first step is to lay out the elements of the approximation table for the transformed 
sequence in the following standard arrangement 



(0) 



E 



E 



(2) 



E 



(3) 



E, 



(0) 



E, 
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(2) 
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So 
Si 

S2 
S3 



E 



(0) 



E 



(1) 



E- 



(2) 



E- 



(3) 



E. 



(0) 



E. 



(1) 



E. 
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E, 



(3) 



E. 



(0) 



E. 



(1) 



E. 



(2) 



E. 



(3) 



(133) 



If we write the exactness requirement of Eq. (|132|) for the quantities s^, . . . , Sn+fc and solve for the 

(n) in) 

unknown s = Ef, , we get for the elements Ef, of the A;th column of the approximation table the 
form 

Sn • • • 'Sn+fc 

gi {n) ... gi{n + k) 



E 



in) 



Qk (n + k) 



1 

91 {n) 
Qk {n) 



1 

gi (n + k) 



Qk {n + k) 



(134) 



Then, making use of standard determinant relations it can be shown that the actual computation 
of the determinants of Eq. (jl34p is not really necessary (this finding is in complete analogy to 
the case of the derivation of the e-algorithm from the Shanks process); in turn, introducing the 
auxiliary quantities 

gi (n) ... gi (n + k) 
91 (n) ... gi (n + k) 



(n) 
9kJ 



9k [n) 



9k {n + k) 



1 

91 {n) 



1 

. . gi{n + k) 



(135) 



9k [n) ... gk {n + k) 

defined themselves in terms of ratios of determinants, a recursive evaluation scheme for both the 
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elements of the table and the g\. ■ can be derived, with the rules 



E 



in) 



n G No, 



(136a) 



'0 




9i{n) , 



N 



(136b) 




,(n) 



E\ 




(136c) 



k-l 



(n+1) _ (n) 3k-l,k 
9k-l,k ~ 9k-l,k 





i G N, i>k+l. 



(136d) 



Using this set of recursive rules, the approximation table can be immediately computed by standard 
techniques once a prescription for the auxiliary sequences gi{n) of Eq. (jl32p is given. 

In fact, many extrapolation schemes can be derived as a special case of the i?-algorithm: among 
them we find (with {xn}^=Q and {j/nj^o being arbitrary sequences which specify some interpola- 
tion scheme): 

1. The Richardson extrapolation process, corresponding to the choice gi{n) = x^, leads to the 
Neville algorithm [see Eq. (3.1.3) of Ref. [83] and Eq. (6.1-5) of Ref . [31]]. 

2. The G transformation corresponds to the choice gi (n) = Xn+i-i [see p. 95 of [25] and 
Ref. [150]]. 

3. The e algorithm corresponds to the choice gi{n) = Asn+i-i- 

4. Certain generalization of the e-algorithm correspond to the choice gi (n) = (s„Ax,i) / Axn, 
being R a difference operator generalizing A and defined by i?'^^^ (s„Ax„) = 
A {R'' (snAxn) /Axn) (see page 108 of [25]). 

5. The p-algorithm corresponds to the choice g2i-i{n) = x~^Sn-, g2i ("-) = a^n*) « = 1, • • • , A:. 

6. The Germain-Bonne transformation [151] entails to the choice gi{n) = (As„)\ 

7. The p-process is given by gi{n) = x„ and gi{n) = As„+j_2 (« G N, i > 2) (see Ref. [18]). 

8. Levin's generalized transformation corresponds to the choice gi (n) = x'^'^A.Sn/yn (see page 



One of the fascinating properties of the E algorithm is the fact that it can be used to accelerate 
the convergence of the sequences for which an asymptotic expansion of the error is known. That is, 
assuming that the error r„ = s„ — s = ai (n) -|- 02 52 {n) + . . . has an asymptotic expansion with 
respect to the asymptotic sequence {gi, g2, ■ ■ ■), i.e. gi+i{n) = o{gi{n)) for i G No, then theorems 
exist (see pp. 68-69 of [25] and [80]) showing that the E algorithm is indeed accelerative for the 
sequence. Moreover, it can be shown for the columns of the extrapolation table (jl33p that -E^"^ 



converges to s faster than El_^: Ej^ — s = o {^E^J^ — sj. The practical importance of this fact 
should not be underestimated, since there are a number of methods known which may help to 
obtain the asymptotic expansions for the error. Again, some of these methods are discussed in 
detail in Refs. [25,80]. 

However, the generality of the approach comes at a price, in view of (i) the difficulty of deriving 
general theorems about the convergence properties of such an acceleration framework and (ii) an 



113 of [25]). 
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enhanced sensitivity to numerical problems, due to the augmented complication of the iterative 
scheme (in fact, some of these problems can be avoided in less general, more specialized and more 
computationally efficient algorithms which may be derived for particular cases). 



2.3 Applications of Sequence Transformations to Convergence Acceleration 
2.3.1 Simple Case: Series arising in plate contact problems 

This review would be incomplete without the discussion of some actual applications of the al- 
gorithms discussed in Sees. 12.2.11 and 12.2.31 The review [31] contains a few very useful practical 
applications with numerical example. Also, we should mention that more numerical examples have 
already been given, e.g., in Refs. [34,148,152]. Here, we discuss a few other examples in depth. In 
the current Section, we start with some easy cases in Sec. 12.3.11 In some sense, the examples dis- 
cussed here complement Refs. [34] and [148] by considering a few cases where a direct application 
of the accelerators is possible in the case of nonalternating series. Specifically, we first consider a 
very simple example which occurs naturally in the context of so-called "plate contact problems" . 
Such problems are well known in mechanics and in engineering applications and concern the con- 
tact between the plates and unilateral supports. While, of course, any discussion on the physics 
of the "plate contact problems" is out of scope of the present report, we just mention that the 
mathematical analysis of these problems may be traced back to the evaluation of special infinite 
series [153], two of which have the following structure: 



The practical use of these series, however, is very embarrassing because of their slow convergence 
for x ~ 1 (for Eq. (jl37p ) or x ~ 6 (for Eq. p38p ). In order to overcome such a slow convergence, a 
number of different techniques have been proposed during the last decade [93, 154-157] . In Ref. [158, 
159], for example, it was shown that the combined nonlinear-condensation transformation (CNCT) 
as discussed in subsection 12.2.61 of this report is able to efficiently accelerate the convergence of 
the series (I137p - (ll38p in problematic parameter regions. In order to illustrate this property, we 
display in Table [2] the calculation of the function Rp{x) for the parameter p = 2 and the argument 
X = 1. For this set of parameters, the series (jl37p are summed by means of the "straightforward" 
term-by-term summation (second column) as well as of the CNCT method (third column). As seen 
from these calculations, the CNCT approach allows us to reproduce the "exact" value of i?2(l) 
with accuracy up to 19 decimal digits for a transformation of order n = 16. In contrast, the partial 
sums i?^"^(l) of 16th order may predict only the two first digits of i?2(l)- In analogy to Table [2] for 
Rp{x), we also present Table [3] for a typical evaluation of Tp{x,b) in a region of slow convergence 
of the series representation ()138p . 

2.3.2 Involved Case: Lerch Transcendent 

As a slightly more involved application of convergence accelerators, we describe the basic strategy 
for the evaluation of the Lerch transcendent [160] 



oo 



,2fc+l 




(137) 




(138) 




x\ <l 




(139) 
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Table 2: Acceleration of the convergence of ii2(l) defined in Eq. (jl37p . In 
the first row, we indicate the quantity n which is the order of the partial 
sum or alternatively the order of the transformation. In the second row, 
the partial sums R^\x) = '^^=0 x^^~^^ /{2k + 1)^ of the input series are 



given while in the third row we indicate the CNCT transforms (jl29p . 



n 


d("-) 


(1) 






1 


1.111 111 111 


111 


111 111 


1.237 516 275 551 033 976 


2 


1.151 111 111 


111 


111 111 


1.233 357 211 766 354 240 


3 


L171 519 274 


376 


417 234 


1.233 701 899 765 Oil 829 


4 


L183 864 953 


388 


762 913 


1.233 700 970 001 046 762 


13 


1.215 850 986 


098 


258 088 


1.233 700 550 136 169 820 


14 


L217 040 046 


740 


350 834 


1.233 700 550 136 169 826 


15 


L218 080 629 


466 


677 578 


1.233 700 550 136 169 827 


16 


L218 998 903 


112 


223 950 


1.233 700 550 136 169 827 


17 


L219 815 229 


642 


836 195 


1.233 700 550 136 169 827 


exact 


1.233 700 550 136 169 827 


1.233 700 550 136 169 827 



This function naturally occurs in the evaluation of statistical distributions, which have application 
in biophysics [148]. We assume x, s, and v to be real rather than complex. 

For negative x, the series (I139P is alternating, and its partial sums can be directly used as input 
data for the delta transformation ([85]) . I.e., for x < 0, we evaluate the transforms 

n 

This approach corresponds to the standard use of nonlinear sequence transformations as efficient 
accelerators for alternating series [31]. 

By contrast, for x > the series (|139|) is nonalternating, and the direct use of the delta transfor- 
mation is not recommended. Additionally, one observes that the series (jl39p is slowly convergent 
for X close to unity. In this parameter region, the use of the CNCT for positive x and x close to 
unity removes the principal numerical difficulties. Thus, for x > close to unity, we evaluate the 
transforms '7cnc('T') as defined in Eq. (|129p . on the basis of the input data given by the partial 
sums Sn as given in Eq. (jl40p . 

For both positive and negative \x\ < 0.3, by contrast, it is computationally advantageous to di- 
rectly sum the terms in the defining power series (jl39p term by term, because of the extremely 
rapid "direct" convergence of the power series in this domain. For negative argument, the series 
is alternating, so one can use the 6 transformation, whereas for positive argument, one uses the 
CNC transformation. In order to illustrate the importance of the correct choice of the acceleration 
algorithm, we display in Table U] the results of the computations of the Lerch transcendent (jl39p 
for the argument x = 0.99 and the parameters s = 1.1 and v = 0.1. As seen from the Table, 
while both, the straightforward term-by-term summation and the 6 algorithm result in a slowly 
convergent result, the 14th order of the CNC transformation allows us to reproduce the exact value 
of $(x = 0.99,5 = 1.1, V = 0.1) with the accuracy of about 10~^^. In contrast, if one evaluates 
the Lerch transcendent for the negative argument x = —0.99 and the same set of parameters s 
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Table 3: Same as Table [2] for r3(l, 1.00000001). 
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1.049 324 


221 


5y5 


503 


955 


209 


376 


363 


455 


1.051 799 867 127 602 653 337 384 365 118 


22 


1.051 681 


688 


751 


729 


032 


310 


253 


771 


189 


1.051 799 780 317 229 791 448 360 992 429 


23 


1.051 691 


320 


524 


361 


301 


753 


922 


155 


807 


1.051 799 780 317 229 791 448 360 994 831 


24 


1.051 699 


820 


379 


948 


685 


582 


513 


817 


375 


1.051 799 780 317 229 791 448 360 995 124 


25 


1.051 707 
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954 


780 


387 


793 


955 


705 


424 


1.051 799 780 317 229 791 448 360 995 160 
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1.051 714 
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259 


950 


1.051 799 780 317 229 791 448 360 995 164 


27 


1.051 720 
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089 


799 
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715 
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1.051 799 780 317 229 791 448 360 995 165 


28 


1.051 725 


486 


189 


637 


836 


692 


682 


642 


537 


1.051 799 780 317 229 791 448 360 995 165 



exact 1.051 799 780 317 229 791 448 360 995 165 1.051 799 780 317 229 791 448 360 995 165 
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Figure 1: Numerical algorithms used for the Lerch transcendent. 

and V, the 5 transformation leads to much faster convergence of the series (I139P if compared to 
the CNCT method [cf. Table [2]. These two simple examples, provides us with a first manifestation 
of the general wisdom that numerical acceleration algorithms have to be adapted to parameter 
domains. 

2.4 Applications of Sequence Transformations to Resummations 
2.4.1 Simple Case: Incomplete Gamma Function 

We follow the same strategy as in Sec. 12.31 As a simple example of a resummation algorithm (see 
Sec. I2.4.ip . we discuss a possible implementation of the incomplete Gamma Function r(0,x) for 
positive argument x, based on the resummation of the divergent large- 2; asymptotic expansion, 
where the series is strictly alternating. Of course, the upshot will be in demonstrating that the 
resummation algorithm permits us to use the divergent large- 2; asymptotic expansion also at small 
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Table 4: Acceleration of the convergence of the Lerch transcendent ^{x, s, v) defined 
in Eq. ()139p . Results are presented for the set of parameters x = 0.99, s = 1.1 and 
f = 0.1. In the first row, we indicate the quantity n which is the order of the partial 
sum or alternatively the order of the transformation. In the second row, the partial 
sums Sji = Y^=Q / + "^Y °f input series are given while in the third and the 
fourth rows we indicate the transformation 5n^{l,SQ) and the CNCT transforms, 
respectively. 



n 


















1 


13.480 


716 


950 


334 


14.323 


956 


240 


313 


16.333 647 069 601 


2 


13.914 


057 


335 


485 


14.872 


118 


432 


283 


16.274 855 206 560 


3 


14.193 


574 


009 


886 


15.328 


655 


340 


349 


16.279 439 130 003 


9 


14.965 


099 


286 


381 


16.183 


960 


654 


046 


16.279 415 474 265 


10 


15.036 


154 


906 


635 


16.214 


534 


170 


121 


16.279 415 474 379 


11 


15.099 


561 


126 


734 


16.235 


397 


170 


154 


16.279 415 474 457 


12 


15.156 


650 


914 


411 


16.249 


604 


756 


607 


16.279 415 474 457 


13 


15.208 


442 


487 


775 


16.259 


256 


744 


711 


16.279 415 474 453 


14 


15.255 


730 


569 


885 


16.265 


802 


155 


733 


16.279 415 474 453 



exact 16.279 415 474 453 16.279 415 474 453 16.279 415 474 453 



positive X with good effect. A much more involved case, which demands an adaptive algorithm, is 
the fully relativistic hydrogen Green function, as discussed in Sec. I2.4.2[ 

We start our discussion on the application of nonlinear sequence transformations to the resumma- 
tion of the divergent series from the computationof the incomplete Gamma function given by [161] 

roc 

r(a,z) = / die"* . (141) 

J z 

For the sake of simplicity, below we restrict our analysis to the special case of a = and of positive, 
real argument: z € R, z > 0. For such z, the incomplete Gamma function reads as 

r(0,z) = j dt — = -Ei{-z) , (142) 

where Ei(z) is the standard exponential integral. 

Usually, the evaluation the incomplete Gamma function (jl42p is traced back to the computation of 
its asymptotic expansion. The explicit form of such an expansion is however different for the cases 
of small and large arguments z. For example, while for large, positive z the incomplete Gamma 
function T(0,z) has the asymptotic expansion 

r(o,^) = ^ ^i-irm - , z^oo, (143) 

^ n=0 ^ ^ ^ 

a completely different series arises for z — > 

2 3 

r(0,z) = lnz-7 + z-^ + ^ + 0(z^), z^O, (144) 
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Table 5: Same as Table H for x = -0.99, s = 1.1 and v = 0.1. 

n Sn (5i°^(l,so) Tc^c{n) 

1 11.697 791 285 548 11.989 386 929"355 9.312 621 614 169 

2 12.131 131 670 698 11.96 5 797 543 541 9.521 513 807 481 

3 11.851 614 996 297 11.967 038 897 324 9.812 306 769 207 



8 12.009 956 315 957 11.967 090 786 6 54 10.403 057 902 502 

9 11.929 460 886 859 11.967 090 786 6 22 10.479 117 052 492 

10 12.000 516 507 113 11.967 090 786 619 10.547 006 982 924 

11 11.937 110 287 014 11.967 090 786 619 10.608 315 040 497 

exact 11.967 090 786 619 11.967 090 786 619 11.967 090 786 619 



where 7 = 0.577216 ... is Euler's constant. It is clear that (11440 has numerical advantages over 
(|143p for small z, whereas for large z, one would certainly refrain from using (|144p . 

As seen from Eqs. (jl43p and (jl44p . asymptotic expansions can be utilized to compute the Gamma 
function r(0, z) in the regions of the large (z — > 00) and small (z 0) positive arguments z. The 
problem arises, however, for the intermediate values of z. In order to resum the divergent series for 
2; — > 00 and, hence, to compute the incomplete Gamma function for such (intermediate) z region 
we may apply the 5 transformation. That is, by using the partial sums of the alternating divergent 
series (|143p 

— Z " / 1 \ ^' 

= VE(-l)'^! U ' (145) 

as input data for Eq. (j85p . it is possible to recover the actual value of the incomplete Gamma 
function to high accuracy, even at relatively small z. In Table El for example, the resummation of 
the asymptotic expansion (jl43p . which is a priori valid only for large z, is performed for the case 
2; = 1. As seen from this Table, while the term-by-term summation of Eq. (jl43p fails to compute 
the Gamma function r(0, z = 1), the 26th-order transformation 6^~^\l, sq) reproduces this values 
up to 16 decimal digits. 



2.4.2 Involved Case: Relativistic Hydrogen Green Function 

In Sec. 12.4. H we have demonstrated how nonlinear sequence transformations may be used for 
the resummation of the divergent series arising from the expansion of the incomplete Gamma 
function. Below another example for the application of the nonlinear sequence transformations 
will be discussed, which concerns the evaluation of the relativistic Green function. This function is 
known to be an essential ingredient in (relativistic) atomic physics and quantum mechanics, because 
it allows us to perform a summation over a complete spectrum of intermediate, virtual atomic 
states. Such a summation is required, for example, for describing quantum electrodynamic radiative 
corrections in atoms (see, e.g., Refs. [35,162-164]), for the computation of hyperpolarizabilities and 
for the evaluation of two- (and multi-)photon processes (see, e.g., Ref. [165]). 

The fully relativistic Green function of an electron bound in a hydrogenlike atom is given in terms 
of the Dirac Hamiltonian i^D = a ■ p + (3 m — Za/\x\, where the a and (3 are standard 4x4- 
matrices [166] in the Dirac representation, p is the momentum operator, Z is the nuclear charge 
number and a is the fine-structure constant. We use units in which h = c = = 1 and write the 
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Table 6: Resummation of the large-z asymptotic series for the function 
r(0,z) = exp(-z)/z Efelo^'/(-^)'' at argument z = 1 [see Eq. (fTi3|) ]. 
In the first row, we indicate the quantity n which is the order of the 
partial sum or alternatively the order of the transformation. In the sec- 
ond row, we list the partial sums s„ = exp(— z)/^ ^^^q A;!/(— z)*' [see 
Eq. (jl45p ]. whereas in the third row, the b transforms are displayed. 
These are defined in ([85]) . 



n 




<5i°^(l,so) 


1 


0.000 000 000 000 000 X 10^ 


0245 252 960 780 961 


2 


0.735 758 882 342 884 x 10° 


0.210 216 823 526 538 


3 


-1.471 517 764 685 769 x 10° 


0220 040 039 579 180 


4 


7.357 588 823 428 846 x 10° 


0.219 508 174 842 629 


22 


3.955 539 479 532 130 x 10^° 


0.219 383 934 395 518 


23 


-9.114 871 523 102 565 x lO^^ 


0.219 383 934 395 519 


24 


2.191 353 397 822 361 x 10^3 


0.219 383 934 395 520 


25 


-5.487 119 942 851 230 x 10^^ 


0.219 383 934 395 520 


26 


1.428 755 174 056 189 x 10^^ 


0.219 383 934 395 520 


exact 


0.219 383 934 395 520 


0.219 383 934 395 520 



Dirac-Coulomb Green function as 



X2 



Xi 



(146) 



Here, z is an energy argument which is assumed to be manifestly complex rather than real values in 
some calculations. A closed- form representation of G{x2,xi, z) is not known. However, it is possible 
and convenient to separate radial and angular variables, and to write the Green function in the 
following 4 X 4-matrix representation [166, 167] for the fully relativistic calculations, 

„/ Gl^{x2,xi,z) Xk{xi)xkHx2) -iGl^{x2,xi,z) Xk{xi)x-k{x2) 
G(x2,xi,z) = 2^ 

\ ^Gf{x2,Xl,z) x'^Ki^l)XK{x2) Gf{x2,Xi,z) X-Ki^l)X-Ki^2) 

(147) 

Here, the radial variables I, and the unit vectors \. The symbol Xk{x) 

denotes the standard two-component Dirac spinor, as given e.g. in Refs. [168,169], 



m 1 Xl 
2 



K + 



2' 



« = (-!) 



j+l+l/2 



j + 



(148) 



and the Cl^_^ 1 are Glebsch— Gordan coefficients in the conventions of Ref. [169]. In Eq. (jl47p . 

the expressions of the form Xk{xi) x'k {X2) are to be understood as dyadic products, i.e. as defining 
2x2 matrices. Of course, k has the meaning of a Dirac angular quantum number. The quantities 
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I and j are the orbital and the total angular momentum quantum numbers, respectively, The 
components of the radial Green function are given in terms of the four quantities G'^{xi,X2, z) 
with a,b = 1,2, and these functions have been historically problematic as far as their numerical 
evaluation is concerned. 

Indeed, as seen from Eq. (|147|) . the evaluation of matrix elements involving the relativistic hydrogen 
Green function entails both angular and radial parts. In many cases, the angular part of the Green 
function can be evaluated using standard, yet sometimes cumbersome angular momentum algebra, 
while the computation of the radial functions G'^^{x2,xi, z) is known to be a very involved task. 
An explicit representation of these four components is given by [see, e.g., Eq. (A. 16) in [170] and 
Eqs. (101), (108) and (109) in [171]]: 

(149a) 
(149b) 
(149c) 
(149d) 



GI^{X2,X1, 




= Q(l + 


z)A. 


_(A,z^, 


2CX2)5_ 


^(A, V, 2cxi 


Gl^{x2,xi, 


^) 


= QcA. 


(A,z/, 


2CX2) 




2CX2) , 


Gf{x2,Xi, 




= QcA+ 


(A,z/, 


2CX2) 




2CX2) , 


Gf{x2,Xi, 


^) 


= Q{i- 


z)A 


|-(A,Z^: 


2CX2)5- 


_ (A, I', 2cx2 



Here, the overall prefactor is 



n- ^ r(A-i/) 

~ 4c2 X2Xi,/x^ r(l + 2A) • 



The quantity c is given in terms of the (complex) energy variable z as c = \/l — z"^, and up to 
order (Za)'^, the quantity A is equivalent to the modulus of the Dirac angular momentum quantum 
number k , 

A = y^?^^^^ = |k| + 0(Za)2 . (151) 

Indeed, it is possible and sometimes even necessary to separate A into an integer part and a 
fractional remainder, 

A = |«|-C, C = rfxT- (152) 

We here identify C as the fractional part of A. Finally, is given in terms of the energy variable z 
as 

^ (Za) . (153) 



z 



2 



The quantities A and B in Eq. p49|) are related to the two solutions, integrable at infinity and at 
the origin, respectively, of the corresponding radial Dirac equation, and given in terms of Whittaker 
M and W functions. The specific form is as follows: 

^±(A, i^, y) = (A - u) M,_^/2,x{y) ± - ^=5) M,+i/2,x{y) , (154a) 

B±i\,i^,y) = (^K + ^^^ ) W,_y2,x{y) ± W,+y2,x{y) • (154b) 

As seen from Eq. (jl49p , the computation of the radial Green functions can be traced back to the 
evaluation of the special Whittaker functions 

M,±i/2,x{y)' W,±i/2,xiy)- (155) 
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For purely real and purely imaginary values of the energy variable z, the quantity c and thus the 
argument of the Whittaker functions remains purely real rather than complex. This situation is 
realized along special integration contours for the complex photon energy [35,163,164,170,172], but 
not universally applicable. We will assume here that in many practical applications, the complex 
argument of y can be made small by choosing a suitable integration contour for the photon energy. 

Here, we thus concentrate on numerical methods, some of which are based on nonlinear sequence 
transformations, which permit the evaluation of the Wcij^^\{x) and Maj^^\{x) functions over a wide 
parameter range relevant for self-energy type calculation [35, 163, 164, 170, 172]. This parameter 
range can be characterized by the first parameter i/ it 1/2 being complex, but small in absolute 
magnitude, and the second parameter A being real, but possibly large in absolute magnitude. 
Pathological cases in which the parameters are close to integers, also have to be treated. 

We rewrite the Whittaker functions W and M in terms of more fundamental hyper geometric U 
and iFi functions The M functions occurring in (jl55p can be written as (see [173], p. 337 and [174], 
p. 296) 



M,±i/2,x{y) = e-^J'y^+^iFi QT^ + A-z.,2A + l,y^ 



= e-^^?/l''l+5-^ iFi(|k| + 1-(5±-C- 1^,2 |k| + 1-2 C,?/) , (156) 

whereas the W functions can be rewritten in terms of the hyper geometric functions U (see [175], 
p. 14) 



W^.±i/2,A(y) = e-^J/y^+^ t/QT^ + A-z.,2A + l,y) 



= e-5S^yl''l+^-^C/(|K| + l-(5± -C-z^,2|k| + l-2C,y) . (157) 
We here define the variable 6± as 

{1 for the case + , 
(158) 
for the case — . 

The transformations (|156p and (jl57p generate iFi and U functions, which are given by 

iFi{a±,bx,y) and U{a±,hx,y). (159) 

The parameters are 

a± = \K\ + l-6±-C-v, 6a = 2|k| + 1 - 2C. (160) 

For the W function, it is possible to rewrite the U function as a 2F0 (see [173], p. 342, [174], p. 317 
and [175], p. 61) 

W^.±i/2,A(y) = e-iyy^^^l\F,{^^\ + \-v,\^\-\-v,-^ 

= e-s?' 2F0 ( |k| + 1 - (5± - C - i^, + 1 - <5± + C - -- 



g-i y .±1/2 (|k| + l-5±-C-v)k{-\^\ + '^-5±+C-v)k .-^g^^ 

k=Q ^ 
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This representation entails the Pochhammer symbol {a)k = r(a + k)/T{a) and constitutes the 
fundamental divergent asymptotic series of the W function for large arguments (i.e., in inverse 
powers of y). Note, however, that because of the factorial growth of the numerators with k, the 
series eventually diverges no matter how large y is. Nevertheless, we anticipate here that it can be 
used as an input for the resummation algorithms introduced in Chap. 12.2.31 We summarize that 
the transformation (jl6ip has given rise to a 2^0 function of arguments and 6^, where 

2-Fo(a±, y) , a± = a± , 6± = + 1 - 5± + C - . (162) 

The parameter b'_^ can assume values very close to negative integers, if the absolute magnitude 
of 1/ and is small. This entails a tremendous loss of numerical significance in high-precision 
calculations, if the asymptotic form of the 2-^0 function ()16ip is used for large \k\ (the problematic 
terms originate for a summation index k = |k| — 1 + 5-1-). In the numerical code, it is therefore crucial 
to separate all parameters of the hypergeometric functions into integer and fractional contributions. 

A 




15 100 250 y 



Figure 2: Numerical evaluation of the W function. 



A last representation useful for small y is the decomposition of the W function in terms of two i-Fi 
functions, 



W,±y2Ay) = e-'^yy^+-^ 

r(2A) 



r(-2A) 



+ 



r(i + i + A 



- ^ 1 - A - z.) 

1 1 

2^2 



1^1 ( ^T^ + ^-^'2A + l,y 



e-3 2/yl'^l + |-C 



r(-2|K| + 2C) 



+ 



T{-\K\ + l-5±-iy + C) 
r{2\K\ -2C) 



2 ' 2 
\ — u,l — 2 X,y 

iFi {\k\ + l-S±-u-C,2\K\ + l-2C,y) 



i-f 1 



k\ + 1-6±-u + C,-2\k\ + 1 + 2C,x) 



(163) 



T{\k\ + 1-6±-u-C) 

One of the two iFi functions in Eq. (I163p has parameters a± and b±, and the other is of the form 
iFi{al,bl,y) , al = -\K\ + l-6±-iy + C, 6± = -2 |k| + 1 + 2 C . (164) 



Again, numerical stability problems can result (both parameters and b'^ can be close to negative 
integers) , and these can also be solved by a separation of the arguments into integer and fractional 
contributions. 

We have thus mapped the problem of the calculation of the radial Dirac-Coulomb wave functions 
given in Eq. (|147p onto the problem of evaluating the M and W functions listed in (jl55p . or 
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of evaluating the iFi and U function given in (jl59p . (jl62p and (jl64p . which is equivalent. The 
parameter range for v required for the calculations reported in Refs. [35, 162-164] is approximately 
given hy u £ (0, n) for the so-called low-energy part (soft virtual photons and by G i (0, Za) for 
the high-energy part (hard virtual photons). The parameter range for v thus is tiny in comparison 
to the parameter ranges for A and y, which both comprise the interval (0, 10^) in the high-precision 
calculations described in Refs. [35,162-164]. 



A 



100 ^ 



Figure 3: Numerical evaluation of the M function. 

For the M and W functions listed in Eq. (jl47p , we can thus devise a numerical algorithm which in 
each case depends only on the value of the second parameter A and of the argument y, not on the 
value of the first parameter v. For the W function, a rather sophisticated differentiated treatment 
of the different parameters is required. All the choices are summarized in Fig. [2j The explanation 
is a follows: 

• Algorithm labeled "i-Fi": We decompose the W function into two hyper geometric i-Fi func- 
tions as in Eq. (jl63p . The two i-Fi functions are evaluated by virtue of their defining power 
series (in the argument y). Suitable truncation criteria for the power series can be constructed 
according to Eqs. (D.3) and D.IO) of Ref. [172] for the first and the second \F\ function on 
the right-hand side of Eq. (jl63p . 

• Algorithms labeled "(jj,'^^" and "5^"^": In an intermediate range of arguments A for the W 
function, we use the divergent asymptotic series on the right-hand side of Eq. (jl6ip as input 
for the sequence transformation (i85l) . 



E 



{\k\ ^\-b±-i-v\ {-\k\ + 1 - (5± + C - 



fc=o ^ ^' 



(165) 



The asymptotic series is resummed either via the nonlinear sequence transformation 

5^(1, so) (algorithm "4°^") (166) 

or 

^f\\,Sn) (algorithm "(5f^") (167) 

where 

n= [Re(|K| -l + (5± -C + z^)l. (168) 

Note that in this notation, k refers to the order of the transformation, which is subsequently 
increased until convergence is reached. By contrast, h refers to the element in the asymptotic 
series beyond which the transformation is started, according to the definition in Eq. (j85p . 
The resummation of the asymptotic series yields values for the W function accurate to at 
least a relative precision of at least 10~^^, in the parameter regions indicated in Fig. [2j 
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• Algorithms "2^*0" and "2-^0": For large y, the asymptotic series (|16ip of the 2-^0 function 
is suitably truncated, according to criteria found in the Eqs. (D.6), (E.14) and (E.15) of 
Ref. [172]. For "2-^0"; the asymptotic series is summed from the first element /c = 1 in 
Eq. (I16ip . For cases where both A and x are large, the asymptotic series is summed in both 
upward and downward directions from the element of maximum absolute magnitude in both 
directions (algorithm "2-^0")- 

The first decisive observation in the context of the current review is that the nonlinear sequence 
transformation (jSSp can be used in order to bridge an intermediate region where a suitable trun- 
cation of the asymptotic series ()16ip fails to produce results of sufficient accuracy, because the 
argument y is not sufficiently large, and yet the power series representation of the W function in 
terms of two i-Fi's according to (I163p also fails because of numerical losses due to cancellations 
among the two confluent hypergeometric functions. The second decisive observation, in the context 
of the current review is that a simple, brute-force implementation of this idea cannot succeed, and 
that considerable further sophistication is required if large parameter ranges are to be covered. In 
particular, the separation into the "J^'^^" and "(5^"^" algorithms, counterintuitive at first glance, con- 
stitutes one of the keys to a numerically effective evaluation of the Dirac-Coulomb Green function 

(HHD. 

For the M function in Eq. (jl55p . we use an evaluation scheme illustrated in Fig. [3l For small argu- 
ments, the M function is evaluated by summing the power series of the iFi function in Eq. (I156p . 
starting with the first element of this power series. This region is denoted as "i-Fi" in Fig. [2j The 
elements of the power series of the confluent hypergeometric function are evaluated recursively in 
analogy to Eq. (D.5) in [172], and a termination criterion analogous to (D.6) in [172] is used. For 
large arguments, the power series is also used, but summation is started with the element of the 
power series with the largest absolute magnitude. This maximum element is determined as outlined 
in Appendix E of [172], and summation in the upward and downward direction is terminated by 
criteria analogous to (E.9) and (E.IO) in [172]. The relevant region is denoted as in Fig. [2j 

In the current application just shown we have discussed the accelerated pointwise evaluation of 
special functions by applying convergence acceleration and summation processes to their defining 
expansions. In the same light, it should be stressed that quite recently [36], the same strategy has 
been proposed for the pointwise evaluation of Fourier series and other expansions (e.g., those into 
orthogonal polynomials). Namely, it has been observed that the application of the simple e algo- 
rithm to Fourier series of functions with discontinuities leads to an acceleration of the convergence 
of the Fourier series, and that the Gibbs "overshooting" phenomenon can be significantly reduced. 
A slight generalization of the problem using complex variables turns out to be beneficial in the 
latter context. One defines the Fourier conjugate series S{t) of a Fourier series S{t) with a zero 
constant term by the following relations. 



S{t) = ^ fflfc cos{k i) + ^ 6fc sin(A; t) , (169) 

k=l k=l 
00 00 

S{t) = ^ fflfe sin(A: t) - ^bk cos{k t) , (170) 



k=l k=l 



F{t)=S{t)+iS{t) = ^(afc-i6fc) exp(iA;t). (171) 



k=l 



In many cases [36], the series S{t) = ReF(t) can be numerically determined to good accuracy by 
applying the e algorithm to the complex series F{t), and selecting the real part of the calculated 
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complex Pade approximants to F{t). Note that the calculations described here for the relativistic 
hydrogen Green function also necessitate the application of convergence accelerators "in the com- 
plex plane," as the energy variable z and, consequently, the parameter A defined in Eq. (jl5ip can 
assume complex values. With these remarks, we would like to close the discussion of the convergence 
accelerators and turn our attention to the summation of divergent series. 

3 Borel Summation 

The contents of the current section on Borel summation can be summarized as follows. 

• In Sec. 13.11 a crude guidance for the reader is provided, and some basic properties of the 
Borel summation process are reviewed. 

• In Sec. 13.31 the performance of a number of resummation algorithms is studied using zero- 
dimensional model theories as example. 

• In Sees. 13. 4| 13.51 and 13.61 we discuss both the mathematical theorems governing the large- 
order behaviour of the perturbation theory of anharmonic oscillators, as well as the practical 
summation of the perturbative expansions, in a number of cases of special interest. 

• In Sec. 13.71 some connections of both summation processes as well as convergence accelerators 
to the renormalization group, and to the determination of critical exponents, is discussed. 

3.1 Quick Start 

Why should there be another review article on Borel summation in quantum mechanics and field 
theory? This question may well be asked in view of a number of excellent review articles on related 
subjects. Here is an incomplete list: The development of large-order perturbation theory with a 
special emphasis on the determination of critical exponents of second-order phase transitions has 
been authoritatively summarized in Refs. [176,177]. The development of generalized perturbative 
expansions with applications to quantum mechanical problems has been described in Ref. [178], 
and further works [179-181]. The generalized perturbative expansions, which are also known as 
resurgent expansions, and their connection to perturbation series in quantum chromodynamics 
has been explored in other works (e.g., Refs. [182-187]). The so-called Lipatov asymptotic form, 
which describes factorially divergent contributions to perturbation series due to the explosion 
of the number of Feynman diagrams in higher orders, has been discussed in Ref. [69]. Finally, 
the significance of so-called renormalons has been elucidated in Ref. [188]: these are additional 
factorially divergent contributions to perturbation series due to a specific type of diagrams whose 
number does not explode in higher orders (but whose numerical coefficients do). 

The answer to the provocative question which was asked in the beginning of the last paragraph is 
as follows: We are trying here to accomplish three goals: 

(i) Our first goal is to establish a connection between the physics literature on the subject, and 
several interesting theorems which have been proven on the mathematical side of the field, 
and to illustrate the concept of distributional Borel summability [189, 190] by a number of 
concrete calculations and mathematical considerations. 

(ii) Our second goal is to describe the current status, and possible further developments, of the 
field in the direction of an enhanced perturbative calculation of critical exponents, where a 
lot of effort has already been placed and where progress eventually relies on a combined effort 
along a careful technical analysis of the problem of parameterizing Feynman amplitudes, on 
the sophisticated use of modern high-performance computers, and on numerical algorithms. 
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(iii) Our third goal is to illustrate various connections of the divergent series to physical processes, 
and to illustrate that a number of phenomena would be less well understood today were it 
not for the fruitful combination of results from mathematics and physics in a field which has 
matured only recently. 

In Appendix \C\ we discuss very briefly how the concept of a resurgent expansion [64] , which 
has already proven to be very useful in quantum mechanics, could be useful to justify certain 
remainder estimates for divergent perturbation series in quantum chromodynamics. Many of the 
sophistications discussed here can be summarized under the unifying concept of the singularities 
of the Borel transform in the "Borel plane" spanned by the argument of the Borel transform. 

Observe that the singularities in the Borel plane have got nothing to do with the so-called Bender- 
Wu cuts which are encountered when the coupling constant of anharmonic problems is continued 
analytically into the complex plane [55-57]. Namely, it has been shown that the energy levels of 
typical anharmonic oscillators (e.g., with a g perturbation) behave as multi-valued functions 
when the coupling constant g is continued analytically into the complex plane: They have branch 
cuts. For the quartic perturbation, the branch points and cuts are known (see Ref. [55]) to occur 
near argg ~ 7r/2 and argg ~ 37r/2. The branch cuts are "short" in the sense that they only extend 
to angles of a few degrees about the quoted complex arguments of g. For a detailed discussion, we 
refer to Refs. [55-57], but a short detour is in order. 

Let us suppose that we start with a particular real g and follow a particular energy eigenvalue 
E = E{g) as a function of complex g, as g \s varied along a path in the complex plane and crosses 
a Bender- Wu branch cut, and let us suppose that g is again varied on the other Riemann sheet 
toward the real axis and to the starting value of g: we end up with a different eigenvalue E(g). In 
some sense, the Bender-Wu branch cuts offer "loopholes" in the complex plane that join the energy 
eigenvalues of the problem. Another curious observation is as follows: suppose we vary g scraping a 
branch cut, directly and exactly through a branch point, perpendicular to the scraped branch cut. 
Let us suppose, furthermore, that we follow the evolution of the functions Ei^ii{g) and Eu^i^g) 
that correspond to crossing from branch cut from sheet I to sheet II and vice versa, across the 
branch point: the surprising conclusion [55] is that the two functions Ei^ii{g) and £'11-^1(5) will 
cross at the branch point, and this constitutes a level crossing as a function of g. These and related 
observations have been crucial for the development of the mathematical theory of the anharmonic 
oscillators and related problems, which is summarized here in the form of a number of theorems. 
In a further context, the investigation of complex coupling constants has inspired the development 
of dilation analyticity [191] (sometimes called "dilatation" analyticity) which has proven to be a 
crucial tool for the calculation of complex resonance energies which include an imaginary part (a 
decay width). Yet, the analytic continuation of g into the complex plane does not always provide 
additional insight in the case of the practical problem of obtaining optimized resummed values 
of perturbation series for a given, real, concrete value of g, based on a finite number of known 
coefficients. Problems of that sort are also investigated in this review. 

In Sec. 12.11 we have given a collection of "recipes" as a guide to the selection of the convergence 
acceleration methods. An analogy of the "recipe" collection for the resummation of divergent series 
is much more difficult to provide, and we will resist the temptation to give one here. The reason 
is that while the discussion in Sec. [2] covers a number of algorithms, we here focus on only one: 
the Borel summation. We trust that the numerous examples treated here will serve as an intuitive 
guidance for the interested reader, when a new problem is encountered. 
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3.2 Selected Methods 



3.2.1 General Formulas 

Since its introduction [45,46], the Borel method has been found to be a very powerful and elegant 
tool in quantum mechanics as well as in quantum field theory and, therefore, has been discussed in 
details in an enormous literature (see, e.g., pp. 233 — 238 of [192], pp. 453 — 456 of [193], pp. 880 — 
887 of [177], pp. 465—466 of [194], pp. 373—376 of [195], pp. 87—92 of [196], pp. 54—56 of [197], 
and pp. 743 — 744 of [198], as well as the seminal work [199]). Let us start our discussion with 
the asymptotic, divergent perturbative expansion of some physical observable f{g) in powers of a 
coupling parameter g. We exchange the argument z of the function defined in Eq. (j22|) by g and 
redefine the perturbative coefficients according to a„ — > a„, as is customary in the literature, 

oo 

J]a„5", 5-0. (172) 

n=0 

Moreover, we assume that the coefficients a„ of this series grow factorially as n — > oo, and that 
f{g) is analytic in a suitable sectorial region of the complex g-plane. Then the Borel transformed 
series 

oo oo 

Bf{u) = y ^^,« = y (173) 

^ ^ n ^ r n + 1) ^ ' 

n=0 n=0 ^ ^ 

converges in some circle around n = by the Cauchy-Hadamard theorem because by assumption 
|an|/n! < A"" for some A > 0. That is, the series ()173p possesses a nonzero radius of convergence 
p > 0, and Bf{u) is analytic in the circle of radius p centered at u = 0. 

We can now make use of the relationship n! = Jq°° u"e~"dti, interchange summation and integration 
so that f{g) can be represented by the following Laplace integral (see for example Theorem 122 
on p. 182 of [6] or p. 141 of [173]), 

fig) = - r eM-u/9)Bf{u)du. (174) 
9 Jo 

The case that the series (jl73p for Bf{u) converges for all u £ [0, oo) actually cannot occur here since 
we explicitly assumed above that the coefficients of the formal power series (j22p for f{g) grow 
factorially as n — > oo. The representation of f{g) by the Laplace integral (jl74p is valid provided 
that Bf{u) admits an analytic continuation to a neighborhood of the positive real semiaxis such 
that Bf{u) grows at most exponentially as n — > oo (see for example p. 141 of [173]), so that the 
Laplace-Borel ()174p integral converges at the upper limit. Then, the Laplace integral (jl74p exists 
for some 5 > and provides a finite expression for f[g) even if the perturbation series (j22p diverges 
for every g ^ 0. 

The Borel summation process can be easily generalized. Let us now assume that a physical quantity 
F{g) admits an asymptotic expansion, 

00 

^(5)~J^an5", 9^0, |a„| < C7a"r(A:n + l), (175) 

n=0 

for some rational number k > 1 and for all integer n with C and a as suitable positive constants. 
Then we define the Borel transform of order k as 

00 

n=0 ^ ' 
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By construction, this series converges in the circle of radius p = a ^ around n = 0. If Bp\u) can 
be analytically continued from its circle of convergence to the whole positive real semiaxis, and the 

1/ k 

continuation is bounded by as u —>■ +00, then F{g) can be expressed by the Laplace integral 
(see [200] or Example 3 on p. 45 and Problem 29(b) on p. 73 of [44]) 



-\ du. (177) 



Under certain conditions, the Borel summation process can even be generalized in such a way that 
perturbation expansions can be summed whose coefficients a„ grow like r(A;n + l) exp(an^/4), with 
k and a being suitable positive constants [201,202]. 

Let us now consider two simple examples, 

00 

A{g) = Y^{-\Tn\g'^ (178) 

n=0 

and 

00 

M{g) = Y,n\g^. (179) 

n=0 

The series A{(j) actually constitutes the so-called Euler series (see Sec. 13.3 of Ref. [31]), which 
has often been used as a paradigmatic example of a factorially divergent, asymptotic series in the 
literature. The Borel transforms have a unit radius of convergence and read 

00 

BAiu) = Y.{-\Tn- = -—, Bm{u) = . (180) 

By inserting the above formulas, which actually constitute analytic continuations for u > 1, into 
the Laplace integral (jl74p we find the Borel sums 

A{g) = i r --P(-^/9)du ^^^^ ^ 1 exp(-u/g)du ^^^^^ 

g Jo i + ^i 9 Jo 1-^ 

The Borel sum of A{g) is well defined and may be expressed as [see, e.g., Eq. (3.02) on p. 40 of [76]] 

A{g) = - eM^/9)Ei{l/g), (182) 
9 

where 

/■oo 

E^{g) = / r^e-'dt, \eLigig)\ < n/2 (183) 

is an exponential integral. 

In the case oi M{g), the analytic continuation (|180p has a pole for g = 1 so that the Laplace-Borel 
integral has a singularity for every g > on the integration path and does not define an analytic 
function of g. In view of these complications, it is not immediately obvious in what sense the 
integral representation (jl8ip can be associated with N{g) and used for its evaluation. One may 
choose a principal- value prescription or encircle the pole in either the positive or the negative sense. 
Indeed, we have [see, e.g., Eq. (3.08) on p. 41 of [76]]: 

Ei{-x±iid) = -Ei(2;) =Fivr, x>0, (184) 

where ^ 

Ei(x) = -PP / r^e^dt (185) 



00 
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is the exponential integral defined as a Cauchy principal value (principal part, PP) for x > 0. The 
notation Ei{—x ziziO) denotes the principal branch Ei{—x) on the upper (+) and lower (— ) side 
of the cut along the negative real axis, respectively. Thus, M{g) defined by the divergent series 
P79p can for g > be expressed as follows: 



We will later see that all these possibilities are actually realized in physically relevant examples. 
The sign of the imaginary part, if it exists, can mostly be fixed on the basis on an immediate 
consideration (e.g., a resonance energy eigenvalue has to describe a decaying, not an exponen- 
tially growing state as time increases). However, there are also situations, as already mentioned 
in Sec. 13.2.11 where the presence of a singularity along the integration path of the Laplace-Borel 
integral indicates the existence of additional contributions that cannot be obtained on the basis of 
perturbation theory alone. 

3.2.2 Mathematical Foundations of Borel Summability 

One of the purposes of this review is to form a bridge between mathematics and physics, by 
discussing exact available mathematical results regarding the Borel summability of specific per- 
turbation series. We will thus need, in the current subsection, concepts recalled in Appendix [Aj 
in particular the notions of the resolvent, the spectrum and the Green function defined in Ap- 
pendix [23] and that of an analytic continuation of the resolvent to the second sheet of its Riemann 
surface, which defines a resonance and is described in Appendix IA.2I 

Let us also recall that for the convergent case, the concept of holomorphy establishes a one-to-one 
correspondence between convergent power series and functions of complex variables: this is the 
Weierstrass notion of analytic functions. If convergent power series are replaced by divergent ones, 
this one-to-one correspondence is in general impossible, because infinitely many functions of a 
complex variable singular at a point (without loss the origin) may admit the same divergent power 
series expansion. Consider the example A{g) given by (jl82p . We know that its series expansion as 
\g\ — > 0, |arg(7| < vr is X^^q^' ("5)"- However, any function A{g) + ^ with q > 0, admits 
the same divergent expansion as \g\ in the same directions. Additional requirements, such 
as uniform remainder estimates, are thus required to set up one-to-one correspondences between 
divergent power series and classes of analytic functions. In any physical problem giving rise to a 
divergent perturbation expansion, the most important thing to check in order to assign the correct 
meaning to the expansion is that this series and the solution fulfill the appropriate additional 
requirements. 

The fundamental result in this connection is a classical theorem of Carleman [203]. Before we state 
the theorem, we need to recall some basic definitions and in particular, the fact that the condition 
Re^"^ > for a complex coupling constant g defines the inner region of a circle of radius R 
centered at Re 5 = R/2, Inig = in the complex (7-plane. 

Definition 1 Let R > and let f{g) be a function analytic in the circle Cr = {^f G C : Reg^^ > 
i?"^} of radius R. We say that a power series '^^=QCLng^ is asymptotic to f , or that f admits a 
'Yl^=o^ng^ o-s an asymptotic expansion (to all orders) as g ^ in Cr if, for each fixed N: 




(principal value prescription) 
. vr (lower contour) 

-ivr (upper contour) 



(186) 





(187) 
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It is easy to check that a function has at most one asymptotic expansion as g ^ in Cr. On 
the contrary, infinitely many functions may admit the same asymptotic expansion in a prescribed 
region. For example, if the function f{g) satisfies Definition [1] above, then any function of the form 
f{g) + Ae~^^^, A £ C, admits as well Yl^=o^n9"' as an asymptotic expansion as g — > in 
because A exp(— has an identically zero asymptotic expansion therein. In order to establish a 
one-to-one correspondence between a function and a (possibly divergent) power series expansion a 
stronger condition than the one given in Definition [1] is needed, which is as follows. 



Definition 2 Let C/j be as above. We say that a function f{g) analytic in C/j obeys a strong 
asymptotic condition and admits Yl'^=o '^n9^ as a strong asymptotic series in Cr if there are A,B> 
such that 



N-l 



fig) - E 



an9 



n=0 



<AB^m\g\^ 



(188) 



We can finally state Carleman's theorem: 



Theorem 3 (Carleman) If YlnLo (^n 9"" is a strong asymptotic series for two analytic functions f 
and g, then f = g. Moreover, there are A and B > such that 

\an\<AB''nl, V n G Nq . (189) 

Proof. See, e.g., Carleman's book [203]. 



The following remarks are in order. Carleman proves also the necessity of the condition (11880 . in 
the sense that no weaker condition can guarantee the result. Functions admitting an asymptotic 
expansion with the above properties are thus determined by the series even if divergent and they are 
therefore called quasi- analytic. Observe that the convergence of a power series yields by definition 
a constructive approximation method for its sum. By the same definition, this cannot be extended 
to divergent series even when the conditions of Carleman's theorem are satisfied. 

Therefore the notion of a regular summation method is introduced: given a pair (/(<?), Yl'^=o 9^) 
fulfilling the conditions of Carleman's theorem, a regular summation method is a constructive 
procedure which makes it possible to recover f{g) out of the divergent power series and yields 
the ordinary sum when the series is convergent. In other words, when applicable, a summation 
method makes it possible to associate a unique analytic function to a given divergent series. In 
general the validity of a regular summation method requires much stronger conditions on f{g) 
and Yln'=o ^n9^ than those of Carleman's theorem. An example is the Stieltjes method, equivalent 
to the convergence of the Fade approximants, which requires very restrictive conditions both on 
/ and on Yln'=o'^"'9^ recalled below. A classical result of Nevanlinna [204], however (see also 
Sokal, Ref. [205]), singles out the Borel method as described above as the one which requires just 
the conditions of Carleman's theorem. We have indeed 



Theorem 4 (Nevanlinna) Under the conditions of Carleman's theorem, the divergent series 
'Yl^=Q'^ng^ is Borel summable to f{g) for g real, < g < R. This means the following: (i) 
The Borel transform Bf{u) = X^n^o ^"""^ ^ priori holomorphic within a circle of convergence of 
radius B^^ , i.e. for \u\ < B^^ where B is given in Eq. Iil89\) . admits an analytic continuation to 
the strip {u G C : |Imu| < T,Reu > 0}. (ii) The Laplace integral 1/g /q°° i?(ii) e""/^ du converges 
absolutely for < g < T and for these values of g one has 

f{g) = - / Bf{u)e~'^'3du. (190) 
9 Jo 
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Proof. For a modern presentation, see [205]. 

In many instances, as we shall explicitly see later on in the examples from quantum physics, it 
is actually possible to verify a stronger statement for Borel summability, valid for functions f{g) 
holomorphic in the sector largg] < ^ + e with e > 0. Namely, we have the following theorem. 

Theorem 5 (Watson-Nevanlinna) Let the pair {f{g), X^^q ^ 5*") f'^^fi^^ following conditions: 
(i) f{g) is holomorphic in the open sector = {5 £ C \ {0} : |arg(7| < | + e}, with e > 0, 

and continuous up to the boundary, (ii) f{g) admits Yl^=Q^n g^ cls a strong asymptotic expansion 
to all orders as \g\ in il7r/2+e) ^-S- there are A > 0, B > such that, for N 

fig) -^ang 

n=0 

Then the series Yl'n'=o ^ng^ is Borel summable to f{g) for all complex g in the closure of the sector 
Of = {ffGC: 0<b| <S-i;|arg5| < e}. 

Proof. This classic result is proven for instance in the book of Hardy [6]. Another proof can be 
found in Ch. XII. 3 of Ref. [44]. 



<AB^Nl\g\^, \g\^0, g e n^/^+e (191) 



As before, all these results can be extended to series diverging as n!" for any a > 1 provided the 
function f{g) is analytic in a wider sector, the rule of thumb being: a sector of opening vr both 
in analyticity and in the validity of the remainder estimate is necessary to absorb any n\. Thus, 
an analyticity sector of opening at least 27r is necessary to absorb a divergence as n!^, and so 
on. Obviously, the analyticity is intended on a Riemann surface sector when the opening angle 
exceeds 27r. We limit ourselves to state the analogue of the Watson-Nevanlinna criterion for the 
Borel summability of order k > 1 , in the sense of Eq. ()177p , remembering that k = 1 corresponding 
to ordinary Borel summability. 



Theorem 6 (Nevanlinna-Leroy) Consider once again the pair {f{g), Xlnio ^9^)> ^^^"^ time under 
the following conditions: (i) For some rational number k > 0, f{g) is analytic in the circle Cr = 
{g G <C : Keg^^/'' > R~^}. (ii) The function f{g) admits S^o'^"^?" ^ strong asymptotic 
expansion to all orders as \g\ — > in Cr, i.e. there are A > 0, B > so that, for iV G N and 
\g\ ^0,5 G Cr: 

N~l 



f{g)-Y^ang'' 



n=0 



< AB^r{kN + l)\g\'^ . 



N 



(192) 



Then the series Yl'n'=o g"" is Borel-Leroy summable of order k to f{g) for g real, < g < R. This 
means that the Borel-Leroy transform of order k. 



Bf\u) 



E 

n=0 



T{kn + 1' 



■ u 



(193) 



a priori holomorphic for \u\ < B ^ , admits an analytic continuation to the strip {ti G C : |Imti| < 
R , Re n > 0} and, for < u < R, one has 



1 



fig) = j^i Bf\u) exp 



9 J 



\ 1/k-l 
9) 



(194) 



where the integral in the right side converges absolutely. 
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As above, in many examples where the perturbation series diverges as fast as {kn)l with A; > 1, 
it is actually possible to verify more, namely the validity of the generalized Watson-Nevanlinna 
criterion, where the analyticity requirement is extended from a disk Cr to an open Riemann surface 
sector. 



Theorem 7 (Watson-Nevanlinna-Leroy) Consider once more the pair if{g),Yln=o^9^'')> ^^^-^ 
time under the following conditions: (i) For some rational number k > 0, f{g) is analytic in 
the open Riemann surface sector = {g ^ CI\{0} • l^'i'gS'l < + T^^J The function f{g) 

admits Yl^=o '^n g^ cls a strong asymptotic expansion to all orders as \g\ — > in ^kTT/2+e! the con- 
dition il92\) holds for \g\ 0, g £ ^kTT/2+e- Then the series Yl'^=o^ng^ is Borel-Leroy summahle 
of order k to f{g) for all complex g in the sector Uf = {5 G C : < |g| < B^^, |arg(?| < e}. This 
means that for < \g\ < , \aigg\ < e one has 



fig) 



— / B) '{u) exp 
kg Jo ^ 



9/ 



l/k' 



\ 1/fc-l 

— I dn . 

9J 
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where the integral in the right side converges absolutely. 



We include the following remarks. For the proof of both theorems, see e.g. [44,205]. Actually, the 
Theorems [6] and [7] follow from Theorems H] and [5] via the change of variable g 1— > g^^^, respectively. 
As already recalled, the Borel summation method is regular, i.e. it yields the ordinary sum when 
applied to a convergent series. Actually it yields more. Given a power series Yl'^=o ^ng" with 
radius of convergence r, the Borel transform Bf{u) = Yl'ri=o ^ entire function, and the 
Borel -Laplace integral 1/g Bf(u)e~'^^^ du converges (uniformly) to the sum f{g) in the Borel 
polygon. The Borel polygon (see e.g. [203] for that notion) always contains the convergence disk, 
but coincides with it only if its boundary is a natural boundary for /, namely if / cannot be 
analytically continued beyond its circle of convergence. In other words, the Borel method applied 
to a convergent series gives the analytic continuation of the sum in a region beyond the circle of 
convergence. This is actually the motivation of its introduction by Emile Borel in 1899 in Ref. [45]. 



3.2.3 Borel-Pade Method 



In Sections l3.2.1l and l3.2.2l we have discussed the general formulas and the mathematical background 
of the Borel resummation approach. Within this approach, the Borel transform (1173P of a power 
series usually has a finite radius of convergence about the origin g = 0. For the evaluation of the 
Borel integral, therefore, Bf[u) has to be continued analytically beyond the radius of convergence. 
Such a continuation, however, requires in principle the explicit knowledge of all coefficients of the 
power series; this situation is rather uncommon in physical applications where only a few series 
coefficients are typically known. Therefore, a Borel method can be a useful numerical tool only if a 
reasonably accurate approximation to the analytic continuation of the Borel transformed series can 
be constructed from a finite number of expansion coefficients. A number of approaches have been 
proposed, therefore, for the approximate analytic continuation of Borel transformed perturbation 
series. For instance, in Ref. [206] it has been argued that the analytical continuation can be achieved 
by evaluating Fade approximants. The first n + 1 terms of the Borel transformed series (|173p can 
be used to construct a diagonal or off-diagonal Fade approximant 

{u) , (196) 



Vn{u) 



In/21 /I(n + 1)/21 
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where the symbol [xl has been defined in Sec. I1.2.1[ We then evaluate the (modified) Borel integral 



TVn{g) = - [ du exp{-u/g) Vn{u) , (197) 
9 JCj 

where the integration contours Cj with j = — 1,0, +1 have been defined in Refs. [207-209]. These 
contours are just straight lines tilted in the complex plane so that the singularities of the Fade 
approximants are appropriately encircled when Jordan's Lemma is applied, to give the pole con- 
tributions that finally yield numerical approximations to the imaginary part of the Borel sum. 
Of course, the result obtained along C_i is the complex conjugate of the results along C+i. The 
arithmetic mean of the results obtained by the integration along these contours is attributed then 
to the contour Cq which involves no imaginary part and can alternatively be obtained by averaging 
over C-i and C+i. The limit of the sequence {TT'ni9)}^=o, if it exists, is assumed to represent the 
complete, physically relevant solution: 

f{g) = hm TPnig) . (198) 

In Refs. [207-209], it has been shown that integrations along the contours Cj have been shown 
to provide the physically correct analytic continuation of resummed perturbation series for those 
cases where the evaluation of the standard Laplace-Borel integral is impossible due to multiple 
branch cuts or due to spurious singularities in view of the finite order of the Fade approximant 
()196p . This is confirmed in the current review via the example of the cubic anharmonic oscillator 
treated in Sec. 13. 5[ 



3.2.4 Conformal Mapping Method 

Apart from the Fade approximants discussed in the previous Section, the analytic continuation 
of the Borel transform can be achieved by various other means. Examples are order dependent 
mappings [210,211] and two-point Fade approximants [212] or Fade-type approximants [213], or 
the first confluent form of the e algorithm [214]. Another possibility consists in a conformal mapping. 
In Ref. [215], it has been found that a combination of a analytic continuation of the Borel plane 
via conformal mapping and a subsequent numerical approximation by Fade approximants leads 
to significant acceleration of the convergence of (slowly convergent) series. Such a combination is 
primarily useful in those cases where only the leading large-order asymptotics of the perturbative 
coefficients a„ in Eq. (I22p are known to sufficient accuracy, and the subleading asymptotics have 
not yet been determined. For such asymptotic series, moreover, we will assume below that the 
coefficients an display an alternating sign pattern. This implies the existence of a single pole 
along the negative real axis corresponding to the leading large-order growth of the perturbative 
coefficients, which we assume for simplicity to be at u = - 1. For Borel transforms which have only 
a single cut in the complex plane which extends from u = -1 to u = —oo, one may employ the 
following conformal mapping: 

z{u) = , (199) 

where z is called the conformal variable. The conformal mapping (jl99p maps the complex n-plane 
with a cut along (—1, — oo) unto the unit circle in the complex z-plane. 

Prom Eq. (|199p we can express the Borel integration variable u as a function z: 

4z 

«(^) = (^31)2 • (200) 
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By expanding now the right-hand side of this relation in powers of z, we can rewrite the mth 
partial sum of the Borel transform (11730 as: 

m 

Bf{u{z)) = Cnz'' + 0(z'"+i) , (201) 

n=0 

where the coefficients C„ as a function of a„ are uniquely determined. Making use of Eq. (j20ip . 
we finally define the (partial sum of the) Borel transform expressed as a function of the conformal 
variable z: 

m 

Bf{z) = Y,Cnz'' , (202) 

n=0 

and apply this transform in order to construct the low-diagonal Fade approximant: 

(z) . (203) 

BJ 

Similar as before, as a last step we have to evaluate the Laplace-Borel integrals along the integration 
contours Cj [see Eq. (jl97p ]: 

T^VM = - [ dueM-u/g)rm{z{u)), (204) 
9 Jc, 

at increasing m and to examine the apparent convergence of the transforms TVm{g)- The appli- 
cation of such a Borel-Pade method combined with a conformal mapping of the Borel plane will 
be discussed in detail in the following Section 13.31 for the calculation of series originating from 
zero-dimensional model theories. 



Vm{z) 



[m/2l/[(m + l)/2l 



3.3 Zero— Dimensional Theory 

Until now we have discussed the basic relations as well as the mathematical background of the 
Borel summation and related methods. Now we are ready to use these methods for the analysis 
of the problems which arise in the different applications of quantum mechanics and quantum field 
theory. In this subsection, for example, we will show how the Borel summation method may help 
us in computing the partition function Z{g) of the zero-dimensional i;A^-theory: 



oo 

1 



-5 $4 



(205) 



where we assume 5 > 0. In fact, this function has been used as a paradigmatic example for the 
divergence of perturbative expansions in quantum field theory (see [177,194]). In particular, it was 
shown that the coefficients in the expansion of this integral in powers of g are equal to the number 
of vacuum polarization diagrams in a (f/^ theory. The series (12050 is an example for the "ordinary" 
Borel summability introduced in Sec. 13.2.21 

By expanding (j205p in powers of the coupling parameter g for 17 — > we may obtain the strictly 
alternating series: 

N=0 N=Q ^ \ ' J 
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Table 7: Evaluation of the perturbation series for Z{g) [zero-dimensional (/)^-tlieory, see 
Eqs. (j205p and (j206|) ] for g = 0.1. In the first column, we indicate the quantity n which 
is the order of the partial sum or alternatively the order of the transformation. In the second 
column, we list the Weniger transforms 6n\i,so), whereas in the third column the results 
of the COMPOS algorithm are displayed (the composition starts with the fourth element 
of the series). Apart from these nonlinear sequence transformations, the Borel-Pade method 
as well as the Borel-Pade method combined with a conformal mapping of the Borel plane 
ZM{g = 0.1) are also used for the resummation of the divergent series. The results of such 
resummations are presented in the fourth and fifth columns, respectively. 



n 


5^(1,^0) 


COMPOS 


TZMig = 0.1) 


T^ZMig = 0.1) 


1 


0890 909 091 




O801 186 280 


0.859 380 049 


2 


0.850 511 945 




O870 363 571 


0865 756 316 


3 


0.856 782 997 




0.851 692 619 


0861 937 975 


4 


0.857 694 858 


0877 198 421 


0.859 335 944 


0.858 100 169 


17 


0.857 608 585 


0.857 608 585 


0.857 608 475 


0.857 608 587 


18 


0.857 608 585 


0.857 608 585 


0.857 608 638 


0.857 608 586 


19 


0.857 608 585 


0.857 608 582 


0.857 608 554 


0.857 608 585 


20 


0.857 608 585 


0.857 608 581 


0.857 608 600 


0.857 608 585 


21 


0.857 608 585 


0.857 608 585 


0.857 608 576 


0.857 608 585 


22 


0.857 608 585 


0.857 608 585 


0.857 608 590 


0.857 608 585 


23 


0.857 608 585 


0.857 608 585 


0.857 608 582 


0.857 608 585 


exact 


0.857 608 585 


0.857 608 585 


0.857 608 585 


0.857 608 585 



By using the known relations for the asymptotics Gamma function for large N [161], it is easy to 
find the large-order asymptotics of the perturbative coefficients Cn- 



Cn 



4^ T{2N + l/2) 
T{N + 1) 



16 



N 



V2- 
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TiN) 
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The factorial growth of the (absolute value of the) coefficients Cn leads to a divergence of the 
perturbative expansion in Eq. p06p . even for small coupling g. Therefore, we have to apply special 
resummation techniques in order to recover the partition function Z{g) from the divergent series 
()207p . In our test calculations, we have applied four resummation techniques: (i) the 5 transfor- 
mation defined in Eq. (j85p . (ii) the COMPOS algorithm [25] discussed in Section 12. H (iii) the 
Borel-Pade method as discussed in Eqs. ()196p - ()198p as well as (iv) the Borel-Pade method com- 
bined with a conformal mapping of the Borel plane [see Eqs. (jl99p - (j204p ]. The results obtained 
from these methods for a coupling parameter ^ = 0.1 are compared in Table [71 As seen from this 
Table, the Borel-Pade approximation is clearly less powerful in the summation of the divergent 
perturbative expansion (j206p than the nonlinear sequence transformations, and for large coupling 
g ~ 10, the Pade method fails. However, if the Borel-Pade method is combined with the conformal 
mapping of the Borel plane, it may significantly accelerate the summation of the divergent series. 
For instance, as seen from the fifth column of the Table[71 such a combined method method appears 
to be as powerful as the nonlinear sequence transformations for the summation of the series (|206p 
taken at the relatively large coupling of g = 0.1. 
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3.4 Even Anharmonic Oscillators 



Mathematical discussion of even anharmonic oscillators. The quantum mechanics of anharmonic 
oscillators in the configuration space is described by the Schrodinger equation ff^V = -^V' i^i 
the Hilbert space L'^{M.^). Here, Hi is the Schrodinger operator whose action (on the appropriate 
domain) is defined by 

1 ^ 1 ^ 

Hi = ^^P^ + ^Y.'^^1^+9 Viiqu ...,qN). (208) 

i=l i=l 

We here assume that the order (. of the expansion about the minimum of the harmonic potential 
is exclusively selected. We obtain by definition an anharmonic oscillator or order in degrees of 
freedom. We can write (I208p as 

Hi = HQ + gVi{qi,...,qN), (209) 

where 

H. = \Y.p" + \Y.-U1 (210) 

i=l 4=1 

is the Schrodinger operator of the A'^-dimensional harmonic oscillator. Here, as usual, 

Pi^-'^^ (211) 
oqi 

is the canonical quantization of the ith momentum, and we still denote by qi the multiplication 
operator quantizing the ith coordinate qi (the reduced Planck constant is h = 1 throughout this 
review). We recall that the eigenvalues of Hq are £'ni,...,r!,jv = ('^i + 1/2) uji + ■ ■ ■ + {njy + 1/2) lon- 

The even anharmonic oscillators in dimension N (or, equivalently, in degrees of freedom) are by 
definition the anharmonic oscillators such that the potential Vi in the classical Hamiltonian (j208p 
is a polynomial in A'^ variables of even degree 2s, s = 2, 3, . . .. The basic property is that, if Vi 
is bounded from below, then it tends to +oo as \\q\\ oo for g > 0. Hence any constant energy 
surface is compact, and any classical motion is localized. Then (see e.g. Theorem Xni.67 of [44]) the 
operator Hi is self-adjoint in the Hilbert space L^(1R"') and its spectrum a{Hi) is discrete. Namely, 
it consists entirely of isolated eigenvalues with finite multiplicity which accumulate at +oo. Since 
we may assume without loss > 0, we can conclude that the eigenvalues of Hi (repeated according 
to their multiplicities) form a positive sequence < EQ{g) < Ei{g) < E2{g) ■ ■ ■ — > +oo, and the 
corresponding eigenvectors form a complete orthonormal basis. The ground state EQ{g) is always 
simple (i.e., non-degenerate). 

The quantity we look at in this context is the standard bound state perturbation theory, or Rayl- 
eigh-Schrodinger perturbation theory, for any eigenvalue Ek[g) around any unperturbed eigenvalue 
Ek{0) = Ek of the harmonic oscillator. The first relevant result is that, no matter how the anhar- 
monic terms are selected to form the perturbation, the Rayleigh-Schrodinger perturbation theory 
always generates a divergent expansion. The reason is that the perturbation is not small with 
respect to the unperturbed operator: no matter how small we can take g, the perturbing potential 
gVg will always overtake the harmonic one for ||g|| large enough. For the sake of simplicity, from now 
on we consider only the one-dimensional case. Here the unperturbed eigenvalues are all simple, and 
thus we avoid the discussion of degenerate perturbation theory. In higher dimension this discussion 
is unavoidable. The results are however the same. For details the reader is referred to [216]. The 
first result is: 



Theorem 8 Let E he any eigenvalue of Hq, and let Vi=2s be an even perturbation. It follows: 
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1. For g > sufficiently small there exists one and only one eigenvalue E{g) of Hi{g) near E, 
limg^oE{g) = E. 

2. The Rayleigh-Schrddinger perturbation theory for E{g) near E exists to all orders, and rep- 
resents a full asymptotic expansion for E{g) as g ^ 0; namely 

M 

E{g)-Y,ang^ = 0{g^'+^), VM G N . (212) 

n=0 

Here, a„ are the coefficients of the perturbations series (n G No/, and oq = E is the unper- 
turbed eigenvalue. 

3. The series is divergent as an ~ [{s — l)n]\ for n — > oo. More precisely, there are constants 
Q < A< B such that 

A"r((s- l)n+ 1) < a„ < B"r((s- l)n + 1) Vn e Nq (213) 



For example, if s = 2, i.e. £ = 4 the series diverges as n!, and if s = 4 it diverges as (3n)!, 
independently of the space dimension. The statement of the theorem summarizes a number of 
contributions. The first proof of the divergence, in one dimension, goes back to Ref. [55]; again in 
one dimension, and for nondegenerate unperturbed eigenvalues the bounds ()213p were proven in 
Ref. [217] and in the general degenerate case in Ref. [216]. 

The question now arises of the summability of the above divergent expansions. There is a first 
sharp distinction between the two following alternative cases: 

1. The eigenvalues of Hq are stable as eigenvalues of the operator family He{g) for g small. 
This means the following: let E be an eigenvalue of Hq, and F any circumference of center E 
and radius r independent of g and smaller than the isolation distance of E (i.e., the distance 
between E and the closest different eigenvalue). Then H(^{g) has exactly one eigenvalue E[g) 
inside F for (7 > suitably small. The stability implies obviously the right continuity of the 
perturbed eigenvalues at (7 = 0: limg^o+ E{g) = E. 

2. The eigenvalues of Hq are unstable as eigenvalues of Hi[g) for g small. In that case, the 
eigenvalues of H(^{g) either disappear into the continuous spectrum, or those inside F are 
at least 2 for any (7 > 0. Notice that instability does not a priori prevent continuity of the 
perturbed eigenvalues as 5 — > 0+. 

In the stable case the Rayleigh-Schrodinger perturbation theory is Borel summable to the eigen- 
values of order s for all finite A^. More precisely: 



Theorem 9 Let Ve=2s be an even perturbation. It follows: 

1. Any eigenvalue E of Hq is stable for g > small enough. 

2. Each eigenvalue E{g) of Hi{g) extends to an analytic function of g in the sector |arg5f| < 
^ ^^2^'*'^ in an s-sheeted Riemann surface. 

3. For any e > 0, there exists some B^{E) so that E{g) and its Rayleigh-Schrddinger perturbation 
expansion fulfill the conditions of the Watson-Nevanlinna criterion (Theoreml^ in any sector 
{5 G C : l^l < B^{E) , [arg^rl < ^^^^^ — e}. As a consequence, by Theorem\^ 
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4- The divergent Rayleigh-Schrd dinger perturbation expansion for E{g) is Borel summahle to 
E{g) for < |g| < B^{E) and |argg| < | - e. 

Proof. Assertion 1 was proven by Simon in Ref. [218], as well as Assertion 2 for = 1 (providing 
a rigorous justification of a discovery of Bender and Wu). Assertion 3 (and hence 4) was proven by 
Graffi and Grecchi in Ref. [219] for the ground state and all A^, and by Hunziker-Pillet [216] for all 
states and all N . 

Here, A^ denotes, as before, the dimension of configuration space. 

Observe that the condition [arg^j < ^ in Assertion 4 effectively ensures that Re 5 > and thus 
the stability of perturbation theory. If A^ = 1, so that Hi reduces to the one-dimensional operator 
Hi=2s = + (f' + 5 9^'^, and much sharper results are available. Namely: 



Theorem 10 Let N = 1 and H2s = p^ + + g q^^ ■ It follows: 

1. For any s > 2 any eigenvalue E{g) is analytic in the cut plane \aTgg\ < vr. 

2. For s = 2, 3, the [N/N] and [N/N — 1] Fade approximants of the Rayleigh-Schrddinger per- 
turbation theory converge to E{g), uniformly on compacts in the cut plane. 

3. For s > 3, the Fade approximants do not converge to E{g). 



Froof. Assertion 1 was proven by Loeffel and Martin [220]. Assertion 2 was proven by Loeffel, 
Martin, Simon and Wightman [221] and represents to the present day the only non-exactly solvable 
example in which the convergence of the Fade approximants has been proven. Assertion 3 was 
proven by Graffi and Grecchi [222]. 

The proof relies on the verification that, for all s, the Stieltjes moment problem (see Sec. I1.2.6P 
generated by the perturbation coefficients {—l)"^ an is solvable, i.e. there is a (positive) measure $ 
on [0, +00) such that 

l-OO 

an= rE"d$(x), 77 G No. (214) 
Jo 

If this is true, all functions 

which are analytic in the cut plane [argf^l < vr, admit the Rayleigh-Schrodinger perturbation 
theory expansion as Taylor expansion at g = (notice that the example M{g) of Section 13.2.11 
is a particular case of (|215p with d<I>(x) = e~^da;). When s = 2,3, the solution <I>(x) of the 
Stieltjes moment problem (j214p is unique. In that case it is known that F^{g) can be expanded in 
a convergent Stieltjes type continued fraction 

F^{9) = (216) 



Oo + 



ai + 



02 + 



03 + ... 



where > and 6^ > for all k. The fraction converges in the cut plane, and the coefficients a/j, 
bk of the continued fraction can be expressed in terms of the coefficients ai. of the series. Now the 
even and odd approximants of the continued fraction (i.e. its even and odd truncations) are just the 
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[N/N — 1] and [A^/A^] Pade approximants to the perturbation series in the sense of Sec. 12.2. 1[ By 
contrast, if s > 3, the moment problem has infinitely many solutions; hence the continued fraction 
diverges and so do the Pade approximants. 

Practical resummations of energy eigenvalues of even anharmonic oscillators. We attempt to verify 
the assertions of Theorems ([8]), ([9]) and (fTOll by concrete calculations, using the quartic anharmonic 
oscillator in the normahzation 

Hig) = Ip' + y+ (217) 
as an example. The Rayleigh-Schrodinger perturbation expansions of its eigenvalues 

^ oo 

EN{g)=N+- + Y,A%g\ (218) 
k=l 

leads to well-known perturbation coefficients A^. E.g., for the ground-state energy EN=Q{g), the 
first six perturbation coefficients are given by: 

- - A"^ - -— A^- — 

^0 - 4 , ^0 - 8 ' ° " 16 ' 

,4 30885 916731 65518401 , , 

Ai = , Af. = , A9, = . (219 

° 128 ° 256 ' ° 1024 ^ ^ 

These coefficients grow factorially and, hence, the perturbation series ()218p for the ground state N 
= of the quartic anharmonic oscillator is divergent. For instance, as seen from the second column 
of the Table El we can verify that a straightforward term-by-term summation of the series ()218p 
at strong coupling g = 1.0 fails to reproduce the exact ground-state energy. In order to find this 
energy, therefore, more powerful resummation techniques have to be used. As it has been proven in 
Theorem I10| for example, one may resum the Rayleigh-Schrodinger perturbation series by means 
of Pade approximants. That is, by making use of the first [(n + 1)/2|| terms of the series (I218p . 

we may construct the Pade approximants ln/2l / |[(n + l)/2]] (g) (we here choose the lower- 

L / Ieo 

diagonal ones, as the higher degree of the denominator polynomial usually leads to a little more 
favourable numerical results) . As seen from the third column of the Table [8l the application of the 
Pade approximant leads to rather slow convergence. For instance, for n = 54, the Pade approximant 
reproduces only four significant digits of the ground state energy of the quartic potential. 

In fact, a much faster convergence can be achieved by combining the Pade approximant with the 
Borel summation method. Since the concept of this combined Borel-Pade method has already been 
discussed in Sec. 13.2.31 here we just recall its basic formulas. That is, the Borel transform of the 
perturbation series (|218p is defined by 

^(^)^o = E^^'' (220) 

k=0 



and later used to construct the Pade approximants 



Vniu) 



In/21 /I(n + 1)/21 



(n) , (221) 

Be,, 



where [xl is defined in Sec. 11.2.11 The resummation is accomplished then by evaluation of the 
(modified) Borel integral along the integration contour Cq introduced in [208]: 

TPn{g) = - [ dueM-u/9)rn{u). (222) 
9 Jcq 
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Table 8: Resummation of the asymptotic series (j218p taken at g = 1.0 for the ground-state energy 
of the quartic anharmonic oscillator (12170 . In the first row, we indicate the quantity n which is the 
order of the perturbation expansion. In the second row, we list the partial sums s„ = Ylk=o ■^n=o9^- 
In the third row, we list the Fade approximants of the perturbation series, whereas in the third 
row, the Borel-Pade transformation is presented [cf. Eqs. (j220p - (j222p ]. 



n 


z_/fc=o o" 


K2l/l(n + l)/2l 


(g) 

En 


TVn{g) 








1 


1.250 000 000 X 10" 


-1.000 000 000 X 10° 


1.657 119 541 


X 


10" 


-1 


2 


-1.375 000 000 X 10° 


6.666 666 667 x 10"^ 


7.197 424 642 


X 


10" 


-1 


3 


1.943 750 000 x 10+^ 


1.542 372 881 


X 10° 


9.643 165 557 


X 


10- 


-1 


4 


-2.218 515 625 x 10+^ 


7.338 453 886 x 


10-1 


7.889 653 148 


X 


10- 


-1 


5 


3.359 128 906 x 10+^ 


9.925 839 156 x 


10-1 


8.122 744 488 


X 


10-^ 


-1 


49 


8.809 956 028 x 10+^^^ 


8.038 179 244 x 


10-1 


8.037 706 514 


X 


10" 


-1 


50 


-1.309 231 884 x 10+^^ 


8.037 397 565 x 


10-1 


8.037 706 511 


X 


10-^ 


-1 


51 


1.984 872 788 x 10+^^^ 


8.038 080 847 x 


10-1 


8.037 706 513 


X 


10- 


-1 


52 


-3.068 686 346 x 10+^^ 


8.037 459 721 x 


10-1 


8.037 706 512 


X 


10- 


-1 


53 


4.836 297 199 x 10+^^ 


8.038 004 311 X 


10-1 


8.037 706 512 


X 


10- 


-1 


54 


-7.767 069 326 x 10+^^ 


8.037 508 557 x 


10-1 


8.037 706 512 


X 


10- 


-1 


exact 


8.037 706 512 x 10"^ 


8.037 706 512 x 10"^ 


8.037 706 512 


X 


10- 


-1 



For this contour, the final result involves no imaginary part. 

By making use of Eqs. (|220p - (|222p . we are now able to carry out a resummation of the pertur- 
bation series (j218p for the ground state of the quartic harmonic oscillator. The results of such a 
resummation, performed for the coupling parameter g( = 1.0 are presented in the third column of 
the Table [8l As seen from this Table, the transform TVn{g) of 53rd order allows us to reproduce 
the exact ground-state energy up to 10 decimal digits demonstrating that with the Borel-Pade 
resummation method one can use the perturbation series at large coupling parameters g. 

As mentioned above, calculations presented in Table [8] have been carried out for the coupling 
parameter g = 1.0. While this parameter already corresponds to the large-coupling regime, one 
often needs to study the energy spectrum of the anharmonic oscillator at even larger couplings. In 
place of a Borel resummation of the perturbation expansion (1218P , for such (very) strong coupling 
g it is usually convenient to perform the so-called Symanzik scaling q i— > g~^^^q in the Eq. ()217p 
and rewrite the quartic potential into another one with the same eigenvalues but a fundamentally 
different structure 



H{g) = g 



1/3 



H^ + ^-^g 



«2 



(223) 



where the large-coupling Hamiltonian i^L does not depend on g: 



H. = -l[£)\,'. (224) 

As seen from these Equations, the eigenvalues of the quartic Hamiltonian for the g ^ oo are 
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ground state (N = 0) 




Coupling parameter g 



Figure 4: Exact values (filled circles) for the ground state energy of the 
quartic anharmonic oscillator as a function of g, together with the partial 
sum of the strong-coupling expansion (solid line) as defined by the first 
three nonvanishing terms in powers of g~'^^l^ [see Eq. ([226])]. 



determined in leading order by the eigenvalues Ej^ of the "unperturbed" Hamiltonian (j224p : 

EM{g) « g^/'^E^^'^ . (225) 

An even more accurate estimate of the eigenvalues Ei\]{g) for large g can be obtained if apart from 
the leading term, higher-order terms are included in asymptotics (|225p . In order to compute these 
higher-order terms we may apply standard Rayleigh-Schrodinger perturbation theory to Eq. (j223p 
and use the fact that the "perturbation" V{g) = q^g~'^/^ /2 is Kato-bounded with respect to the 
unperturbed Hamiltonian H]^ for large g. Within such an approach, the convergent strong-coupling 
perturbation expansion can be written for each energy E]<s(g): 

oo 

EN{g)=g'/'Y.^N9-'''^'. (226) 

fc=0 

Here, we have = E^^\ and the higher perturbation coefficients L^^^ can be calculated to high 
numerical accuracy in the basis of the wavefunctions of the unperturbed Hamiltonian ()224p . As 
seen from Fig. [H the asymptotic expansion (j226p as defined by first three nonvanishing terms Lj], 
Lq and Lq is able to reproduce with a good accuracy the exact energy of the ground state of quartic 
Hamiltonian (j217p starting with relatively moderate coupling oi g = 0.4. 

3.5 Odd Anharmonic Oscillators 

Mathematical Discussion of odd anharmonic oscillators. The odd anharmonic oscillators in dimen- 
sion N (or, equivalently, in N degrees of freedom) are by definition the anharmonic oscillators 
such that the potential Vg, in the Hamiltonian (|2U8p is a polynomial in variables of odd degree 
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2s + 1, s = 1, 2, . . .. The most famous example is probably the potential V = q2, yielding the 
two-dimensional Henon-Heiles Hamiltonian 



n=^{pl+ pI) +ujfql+ ujI ql + g qj , (227) 

which represents one of the fundamental models for the onset of chaos in classical Hamiltonian 
mechanics. 

The basic difference to the even case is that here constant energy surfaces are in general not 
compact. Then the classical motions reach infinity when the initial conditions belong to open 
regions in phase space, and actually do it in a finite time. The quantum aspect of this mechanical 
property is the lack of self-adjointness of the corresponding Schrodinger operator, and consequently 
the lack of uniqueness of quantum dynamics. To better understand this point, consider the simplest 
one-dimensional case, the cubic anharmonic oscillator, whose Hamiltonian is (in an appropriate 
normalization) : 

Hi=3=P^ + q^ + 9q^- (228) 

The potential Vi=^{q) = + gq^ has a minimum at g = with value V3(0) = and a maximum at 
q = —Ij'ig with value V3(— 2/3(7) = Therefore, all classical trajectories with initial conditions 

fulfilling E > 2^ reach — cxd in a finite time (with at most one inversion of velocity) because the 
time required is well known to be 

(229) 



sjE - cp' - gq 



3 



and this integral is convergent. Under these conditions (the classically incomplete case: see e.g. [44]) 
the corresponding Schrodinger operator ip^^q^^gt^ is not essentially self-adjoint because there are 
infinitely many different ways in which it can be realized as a self-adjoint operator in the Hilbert 
space L^(IR). On the other hand the Rayleigh-Schrodinger perturbation theory exists. Actually the 
perturbation theory of the cubic anharmonic oscillator represents the very first problem considered 
by Heisenberg with matrix mechanics: see e.g. Ref. [223], and the first obvious question is to isolate 
the quantity to associate with it. 

We consider this problem only in one dimension, and for a general interaction of the type g^*"*"^, 
i.e. we consider the Schrodinger operators 

i/2s+i(<7) + +59''+', 5 = 1,2,3..., (230) 

The lack of self-adjointness for 5 7^ entails a fortiori the instability of the harmonic oscillator 
eigenvalues. The asymptotic power series E{(j) ~ ^Y^=Q^ri9^ generated by Rayleigh-Schrodinger 
perturbation theory near any unperturbed eigenvalue E = ao has indeed the following properties: 

1. It contains only even powers; namely, 02^+1 = 0. This is because the unperturbed operator 
is even under the parity operation {Pip){q) = ip{—q), and the perturbation is odd. 

2. The nonzero coefficients a2n have constant sign, the perturbation series is nonalternating. 

3. The even coefficients behave as a2n ~ (sn)] for n — > 00. 

The proof is given in Ref. [224]. If the series has only even powers of constant sign, the series 
^'^=0 CLn 9)^ has alternating signs. It is thus conceivable that, in examining the summability 
properties of Rayleigh-Schrodinger perturbation theory, the natural starting point should be the 
examination of the operator H2s+i{g) for g complex instead of g real. One has indeed 
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Theorem 11 Let H2s+i{g) he the Hamiltonian of an odd anharmonic oscillator as discussed. It 
follows: 

1. If |arg/3| < ^, the operator i?2s+i(i/9) can he realized in a unique way as a closed operator 
in L^(]R) with discrete spectrum; moreover (i?^2s+i)*(i/3) = -f^2s+i(— where f3 denotes the 
complex conjugate of (3. 

2. The eigenvalues of the unperturhed operator Hq are stahle as eigenvalues E{\(3) of the argu- 
ment i/3 for \(5\ small and |arg/3| < ^; that is, if Eq is any eigenvalue of the one- dimensional 
harmonic oscillator (and hence simple) there is B{Eq) > such that II{\(3) has exactly one 
eigenvalue E{\f5) near Eq for \ f3\ < B{Eq) and |arg/?| < |. 

3. Any eigenvalue E{iP) is a holomorphic function of (3 in the region {(3 : < < 
B{Eq); |arg/3| < the function E{\f3) admits an analytic continuation to any sector 
Ss = {13 : Q < \(3\ < BsiEo) ; -ZLiH£zll + ^ < arg/3 < 1^2^ - 6} , \f 6 > on a {2s + 3)- 
sheeted Riemann surface. 

4. The function E{\0) and the Rayleigh-Schrddinger perturhation series fulfill the conditions of 
the Watson-Nevanlinna Theorem [5| in any sector {(3 : Q <\(3\ < Bs{Eo) , -ZLiHpii + 5 < 
arg/3 < ^(H|tl) v<5 > 0. 

5. The Rayleigh-Schrddinger perturhation series is Borel summahle to E{\(3) in any sector {(3 : 
< < Bs{Eq), |arg/3| < f — 5}, ^6 > or equivalently, to E{g) in any sector {\g\ < 
Bs{Eo) ; 6 < arg g < ir — 6}, if we set g = \(3. 

Proof. These results are proven in [189,190]. 

The above results clarify the meaning of the perturbation series as long as the summation is 
performed along nonreal directions; if the Borel sum is continued to g real the result is a complex 
function. The question thus arises of establishing the relation between this complex-valued function 
and the real-valued perturbation series, and of finding the meaning of the imaginary part. The 
physical intuition is that this complex valued function should represent a resonance of the problem, 
even though there is no continuous spectrum involved: the real part should be the position, and 
the imaginary part the width, because it should be proportional to the penetration of the barrier. 
The distributional Borel summability puts on a rigorous basis this intuition. Namely we have: 

Theorem 12 Let E{g) he the analytic continuation to Img = of the function E{g), which exists 
hy the ahove Theorem. Then: 

1. E{g) and the perturhation series X]^o^'"fi'" fulfill the conditions of the distributional 
Nevanlinna-Leroy criterion (see Theorem\2^ and Definition\2^ in Appendix B) of order ^"^^ . 

2. The divergent perturbation series Yl'n'=o 5" is Borel summahle in the distributional sense 
of order s to E[g) for \g\ suitably small. Moreover, IieE{g) is even in g and Im E{g) is odd. 

Proof See Ref. [225]. 

A few remarks are in order, (i) Equation (1.6) of Ref. [225] clarifies the role of the standard WKB 
analysis to the barrier penetration and the imaginary part Im E{g), in the context of odd anhar- 
monic oscillators, (ii) The analogue of Theorem [12] for the Henon-Heiles Hamiltonian (3.4.1) and 
for the ground state energy level of the harmonic oscillator has been recently proven in Ref. [226]. 
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(iii) Consider the odd oscillators for purely imaginary values of the coupling g, i.e. the operators 
H{i P), (3 G R. Let P be the standard parity operation: {Pip){x) = ip{—x)^ and T the time-inversion 
operation, equivalent, as we know, to complex conjugation: {Ti/j){x) = ip{x). P and T are unitary 
operations, and so is their product PT = TP. Now it is immediate to check that 

[H{i(5),PT] = [H{ip),TP] =0, (231) 

i.e. the odd anharmonic oscillators are examples of PT-symmetric operators. It has been con- 
jectured around 1985 by Bessis and Zinn-Justin that the spectrum of the odd PT-symmetric 
anharmonic oscillators should consist of real eigenvalues for all values of /? £ IR, even if the op- 
erator H{(3) is obviously not self-adjoint and not even normal because it is easily verified that 
[H{if3),H{il3)*] = [H{if3),H{-if3)] / 0. Note that by the Borel summability the reality of the 
spectrum is valid for (3 small (depending on the eigenvalue index k) because the eigenvalues of the 
harmonic oscillators are stable for /3 small and the perturbation series for each eigenvalue is real 
when g = i /3, /3 G R. Notice that in this general setting, the non-normality makes also the prob- 
lem of the existence of eigenvalues a nontrivial problem in situations where perturbation theory 
cannot be applied, such as the one obtained for /3 large. Later Bender (Ref. [227]) conjectured that 
PT-symmetric Hamiltonians should always have real spectrum, so that PT-symmetry could be con- 
sidered a condition equivalent to self-adjointness in ensuring the reality of the spectrum as long as 
Schrodinger operators are concerned. This latter conjecture has been disproven [228], but the origi- 
nal Bessis-Zinn- Justin conjecture has been recently proven by Shin [229] (see also Refs. [230-232]). 

Theorem 13 Let H2s+i{g) be the Hamiltonian of the one- dimensional odd anharmonic oscillator. 
For any s G N and any (7 = i /3, /3 G R, the spectrum of the operator H2s+i{g) consists of countably 
many simple positive eigenvalues Ei^{\(3). 

Proof. See Shin [229]. 

Finally, we note that similar results as for anharmonic oscillators hold for the resonances of the 
Stark effect, see Appendix IB. 2[ 

Practical resummations and complex scaling of eigenvalues for the cubic potential. As for the even 
anharmonic oscillators, we now supplement the mathematically inspired analysis by some concrete 
calculations. 

Before starting the numerical analysis of the cubic potential, however, we may mention that the 
problem of "recovering" the complex resonance energies from a nonalternating perturbation series 
has been discussed for the Stark effect in [207, 209] and for the quantum electrodynamics (QED) 
effective action in [208]. For the Stark effect, for example, the situation is as follows: the problem 
is spectrally unstable: as soon as the electric field is turned on, all eigenvalues of the hydrogen 
atom disappear in the continuous spectrum (a phenomenon first proven by Titchmarsh [233]). 
Thus the perturbation theory is always divergent. The physical reason is that the linear potential 
of the electric field generates a tunneling phenomenon. As first pointed out by Oppenheimer [234], 
the hydrogen eigenvalues are thus expected to turn into resonances: their width (inverse life-time) 
should be proportional to the probability of penetrating the barrier. In the 1970s and 1980s, this 
picture has been set on rigorous mathematical grounds: moreover, the perturbation theory of the 
problem has fully been understood. Not only it can be uniquely associated with the resonances via 
ordinary and distributional Borel summability, but it defines very effective approximation schemes 
for computing both the location and the width of the resonances. 

Here, we can rely on the literature for the Stark effect and for the QED effective action and 
treat the conceptually simpler case of the cubic anharmonic oscillator as an example. Calculations 



66 



for such a potential will be performed in order to illustrate how to derive the complex resonance 
eigenvalues (including the width) from a real perturbation series using Borel summation in complex 
directions of the parameters, in the sense of distributional Borel summability. We discuss the 
problem of determination of the complex resonance energies for the case of the cubic oscillator, in 
the normalization: 

H{g)=^ + ^ + ^q'. (232) 

The energy spectrum of the anharmonic oscillator can be analyzed numerically, e.g. via a diago- 
nalization of the Hermitian matrix within a (large) basis of harmonic oscillator wavefunctions. One 
has to take into account, though, that the bound states of the (odd) harmonic oscillator become 
resonances, labeled by an index N according to the above Theorems 1111 and 1121 

EN{g) = ReEM{g)-i^^^, (233) 

where T]\f(g) is the width. An appropriate ansatz for studying the energy spectra of the odd 
anharmonic oscillators has been proposed in Refs. [191] and [235] and is based on the complex (or 
so-called dilation) transformation of the coordinate q: 

q^qe'\ (234) 

where 6 may be real or complex. By inserting this transformed coordinate into Eq. (j232p . we may 
obtain the complex-scaled Hamiltonian for the cubic oscillator: 



-2ie I ^„2i0 
2dq^ 2 



e-^^^(-i^ + ^e^^^ + ^gV^^) . (235) 



As discussed in Ref. [235], by finding the eigenvalues of this dilationally transformed Hamiltonian 
one may find that when 9 becomes large enough, one of the rotated branches passes the position 
of a resonance pole, and a new eigenvalue of the Hamiltonian (j235p appears at the position of the 
resonance pole. This resonance eigenvalue, once uncovered, remains fixed in the complex plane and 
is independent of further increases in 9. Hence, the evaluation of the eigenvalues of the Hamiltonian 
(|235p allows us to determine the (complex) energies of the original cubic Hamiltonian (j232p . 

Now we are ready to calculate the energy values of the complex-transformed (I235p and, hence, 
of the original cubic Hamiltonian (j232p . Again, these calculations can be performed via matrix 
diagonalization in the basis of harmonic oscillator wavefunctions. The result of such diagonalization 
procedure for a coupling parameter g = 0.01 is presented in the last line of the Tableland is found 
in perfect agreement with previous calculations [236]. 

We not study the resummation of the perturbation series within the complex plane. In order 
to discuss this resummation procedure let us again recall that within the standard Rayleigh- 
Schrodinger perturbation approach, the eigenvalues of the cubic potential are given by the series 

^ oo 

EN = N + - + Y,A%g\ (236) 
fc=i 

where the first six perturbation coefficients for the ground state = are: 
Ai - 11 42 _ 456 3 _ 39709 

^4 _ 19250805 ^5 _ 2944491879 _ 1075012067865 

° ~ 2048 ' °~ 8192 ' °~ 65536 ' ^ ^ 
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Table 9: Resummation of the asymptotic series (|218p taken at g = 0.01 for the ground-state energy 
of the cubic anharmonic oscillator (I232p . In the first row, we indicate the quantity n which is the 
order of the perturbation expansion. In the second row, we list the partial sums Sn = Ylk=o ^off'^ 
of the purely real perturbation theory, whereas in the fourth row, the Borel-Pade transformation 
is presented [cf. Eqs. (I220]l -(|222D]. 



n 






i 


0.486 250 


000 


0.486 949 910 - 


0.000 000 000 i 


2 


0.484 796 


875 


0.484 497 987 - 


0.000 000 093 i 


3 


0.484 486 


648 


0.484 341 583 - 


0.000 004 008 i 


4 


0.484 392 


650 


0.484 315 106 - 


0.000 010 069 i 


5 


0.484 356 


707 


0.484 315 279 - 


0.000 009 966 i 


10 


0.484 320 


462 


0.484 316 006 - 


0.000 008 057 i 


11 


0.484 318 


395 


0.484 316 001 - 


0.000 008 052 i 


12 


0.484 316 


573 


0.484 315 999 - 


0.000 008 058 i 


13 


0.484 314 


836 


0.484 315 998 - 


0.000 008 058 i 


14 


0.484 313 


053 


0.484 315 999 - 


0.000 008 059 i 


15 


0.484 311 


093 


0.484 315 997 - 


0.000 008 060 i 


16 


0.484 308 


792 


0.484 315 997 - 


0.000 008 060 i 


exact 


0.484 315 997 


0.484 315 997 - 


0.000 008 060 i 



Similar to the case of the quartic oscillator, this series is divergent and, hence, its term-by-term 
summation cannot reproduce the energy values of the cubic potential. Apart from such a non- 
summability, moreover, another problem should be mentioned. Namely, no complex part (i.e. width) 
of the eigenvalues can be obtained by straightforward summation of Eq. ()236p . We may calculate 
these widths, however, if we apply the Borel-Pade resummation method with an integration contour 
deformed into the complex plane as discussed in Eqs. ()220p - ()222p . For such an integration, however, 
we have to choose the proper contour in the complex plane in order to reproduce the correct sign 
for the imaginary part of the energy. In our calculations for the cubic potential, for example, the 
contour C_i has been chosen because it leads to the correct negative sign of the imaginary part 
of the resonance energy eigenvalue. Notice that the distributional Borel summability also does not 
fix the sign of imaginary part on the basis of the purely real perturbation series (see Theorem [T2]) , 
and the choice of the contour has to be inserted based on an "external," physical consideration. 
As seen from the third column of the Tabled the Borel-Pade method allows us to reproduce the 
complex resonance energy of the ground state of cubic potential for a coupling parameter g = 0.01. 

3.6 Double-Well and Fokker-Planck Potential 

Mathematical discussion of double-well like potentials. Let us restrict ourselves to the discussion 
of the one-dimensional case. The simplest example is represented by the symmetric double-well 
Hamiltonian, which reads (in an appropriate normalization) 

H^^{g)=p^ + q\l-gqf (238) 
The translation q ^ q + 1/2^ maps the Hamiltonian to the form + S'^((7 + l/25)^((? — l/2g)'^ 
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which is manifestly an even function with two symmetric quadratic minima at q = ±1/2(7. The 
Hamiltonian (I238|) has minima at q = and q = 1/g and the second one disappears as (7 — > 0. 

Within even anharmonic oscillators examples exhibiting spectral instability have a special status 
and a typical examples is just the double-well oscillator. As already mentioned, the basic phe- 
nomenon in this context is the spectral instability, or, equivalently, asymptotic degeneracy of the 
eigenvalues. For the symmetric double well, the well known physical picture goes as follows (see 
e.g. [237]): in classical mechanics the two wells are identical, and therefore generate the same mo- 
tions, in the sense that the phase curves in the two separate wells are the same. Quantum mechanics 
just adds to this picture the tunneling phenomenon. Therefore we may expect to find, for g small, 
a doublet of eigenvalues of the double well operator {g) close to each unperturbed eigenvalue 
E, and therefore E is unstable. The first mathematical proof of this picture is due to Harrell [238] 
who showed that as [^l — > 0, there are exactly two eigenvalues E±{g) of H^^{g) near any unper- 
turbed eigenvalue E, and that the Rayleigh-Schrodinger perturbation expansion is an asymptotic 
expansion for both eigenvalues of the doublet. 

As we have seen, to understand the meaning of perturbation theory, it is in general necessary to 
understand the analyticity property of the eigenvalues for complex values of the coupling constant 
g. In the case of only a single well, this is possible thanks to the stability of the eigenvalues even 
for complex values of the coupling constant. It is remarkable that also in the unstable case of the 
double- well, one can identify and construct a pair of operators H±{g), the so-called single well 
operators which admit E±{g) as stable eigenvalues as g — > in a suitable complex region. The 
operators H^{g) and H-{g) are obtained by projecting orthogonally H{g) onto the subspace of 
the even and odd functions with respect to the top of the barrier at x = l/2g, respectively. This 
represents the mathematical justification of the physical intuition that the two perturbed quantum 
states "live" in two separate but equal wells, and makes it possible to perform the stability analysis 
as already mentioned. The result is 

Theorem 14 Consider the double well operator H^^{g). 

1. 3R > so that both eigenvalues E±{g) are analytic in the Nevanlinna disk Cr = {g : 



2. The eigenvalues E±{g) are stable near the eigenvalue E of the unperturbed problem as \g\ — > 0, 
g G C/j, if considered as eigenvalues of the operators H±{g), respectively. 

Proof. See [190,239,240]. 

The summability statement for the perturbation theory can hold only in an indirect way, identified 
in [240]. Two further operators have to be associated with the differential expression p^+g^(l — (7 g)^, 
namely the resonance operators T±i{g). The operator Tj^i{g) is defined by the the differential 
expression + q^{l — qq)"^ by modifying the boundary condition at +00: the functions in its 
domain are required to be entire analytic and vanish as \q\ +00 within the sector 5i = {g G C : 
^ < argg < f}- Analogously, the operator T^i{g) is defined by the vanishing boundary condition 
in the symmetrical sector S'_i = {xGC:— ^< argj; < —^}. These considerations allow us to 
express the matrix elements of the resolvent of (g) as linear combinations of those of T±i. 

We now define various projection operators. Let 



represent the projections onto the subspaces of the even (respectively odd) functions with respect 
to the barrier point x = l/2g. Then, H^'^{g) = P±{g) {g) P±{g). We also define the following 





(239) 
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projectors, 

= 2,i-£_^^^^-^d.. (240b) 

These project out those eigenfunctions oi (g) and H^'^{g) whose energies fulfih \En — E\ < 
r. Let ipo be the unperturbed eigenvector corresponding to E. With the convention tjj± = P±{g) V'Oj 
we now define the expressions 

N±ig)= {i^o\HY^i9)Q±i9)\i^o) = {i'±\H'''^{9)Qi9)\^±) , (241a) 

D±{g) = (^0 \Q±{9)\ ^o) = . (241b) 

By standard Rayleigh-Schrodinger perturbation theory, we can now write any double well eigen- 
value E± (g) under the form 

"""^^^ -D^)- (^.IQ(.)I^.) • ^'"'^ 

If we restrict our attention now to a specific manifold of states with a definite parity, then we may 
write the energy eigenvalue of the symmetric double-well as 

Ete) = ||. (243) 

Following Ref. [240], it is possible to prove a representation formula of the type 

N{9) = lM9) + m] (244) 
and analogously for D{g), g £ Cr. Then one proves: 



Theorem 15 With the above definitions, we have the assertions: 

1. The function ^{g) admits a power series expansion in g; ^{g) and its expansion fulfill the 
conditions of Theorem\2^ (see Appendix B). More precisely: 

2. The power series expansion is Borel summahle in the distributional sense (see Definition \20[ 
Theorem\2^ and the remark after it) i.e.: ^{g) is the upper sum, ^{g) is the lower sum, N{g) 
is the distributional Borel sum, and ^{g) — ^{g) is the discontinuity. 

3. An analogous statement holds for D{g) and its perturbation expansion. 



Proof. See [240]. 

Let us add a couple of remarks concerning the effect that the asymptotic degeneracy phenomenon 
and the related tunneling are fragile, i.e. they critically depend on the exact symmetry. Indeed 
these phenomena disappear when the the symmetry is broken even by an infinitesimal amount. 
The following two observations can be made, (i) Consider for example the Schrodinger operator 
Hg ,: = + — 2g{l + e)q^ + g^ q'^, e > which corresponds to shifting the second minimum up 
by e from 0. Then the unperturbed eigenvalues are stable and the perturbation expansion near any 
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unperturbed eigenvalue is Borel summable (for the proof, see [241]). (ii) Suppose we perform an 
arbitrarily small variation of the profile of the potential g^(l — gqY away from the minima. Then 
the constant A changes by a finite amount, so that the tunneling probability changes by several 
orders of magnitudes (see for instance [217,242]). The latter observation has been termed the "flea 
over the elephant" by Simon in Ref. [243]. 

A remarkable example by Herbst and Simon [244] shows that this phenomenon takes place even 
within perturbation theory of stable, but asymmetric multi-well anharmonic oscillators. Consider 
indeed the following Schrodinger operator 



where the symmetry of the minima is broken by the linear term. This symmetry breaking makes 
the unperturbed eigenvalues stable; however the Rayleigh-Schrodinger perturbation theory for the 
ground state E{g) is identically zero, while E[g) itself is not zero. 

This problem finds a natural solution here within the context of resurgent expansions and has also 
been treated in [240], where a general Schrodinger operator is examined of the type 



which reduces to the symmetrical double well for j = and to the Fokker-Planck Hamilto- 
nian [181, 244] for j = —2. The precise mathematical statement of the solution is too technical 
to be reproduced here. Roughly speaking, however, the results are as follows: 

(i) If j 7^ 2p, p = 0, ±1, . . ., the same results as Theorem 1151 hold . More precisely, the eigenvalues of 
H{g,j) split into two disjoint subsets {Ej\f{g,j)} and {E'j^{g, j)}. The eigenvalues {E]\f{g,j)} are 
stable near the eigenvalues 2N + 1 + j'/2 of H{g,j), the harmonic oscillator up to the shift j'/2; 
The eigenvalues {E'j^{g, j)} converge to 2N + 1 — j/2 and therefore have nothing to do with the 
limiting problem. They are instead stable with respect to another limiting problem at g = 0, namely 
the harmonic oscillator up to the shift —j/2. The eigenvalues Ei^{g,j) admit a Borel summable 
perturbation expansion, while the perturbation expansion for any single eigenvalue E'j^{g,j) is 
identically zero, (ii) If j = 2p, p = ibl,ib2... an eigenvalue of the first limiting problem may 
coincide with an eigenvalue of the second one. Therefore a 2 x 2 matrix has to be introduced, with 
elements still defined in terms of matrix elements of the resolvents of the corresponding single well 
operators. The two perturbed eigenvalues are the eigenvalues of this matrix. The coefficients of the 
matrix have expression analogous to (j244p . and are likewise Borel summable in the distributional 
sense. 

Practical resummations of double-well like potentials. We start from the case of the symmetric 
double-well potential, this time in the normalization: 



where, again, 5 is a real, positive coupling parameter. Obviously, if g = 0, then Eq. (j247p represents 
the Hamiltonian of the quantum harmonic oscillator whose eigenvalues are given by the well known 
formula E'tv = N + 1/2, where G Nq is the principal quantum number. For nonvanishing 
coupling, by contrast, no analytic expression for the energies Eiy{g ^ 0) can be derived and, hence, 
approximate methods have to be applied for solving Eq. (j247p . First, we again start from the 
Rayleigh-Schrodinger perturbation theory which allows us to write the perturbation series (j218p 
series for the energy values E]y{g). For the ground state N = 0, the first six coefficients ^^=0 
these series are given by (see, e.g., [85]): 



Hppig) =p' + q\l-gqf + 2gq-l 



(245) 



H{g,i) = / + q\l - gqf - j{gq - 1/2) 



(246) 



(247) 



- -1 42 _ _Z /l3 _ 

^0 — ' ^0 — 2 ' — 



89 



2 




5013 



88251 



46 



3662169 



(248) 
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Table 10: Energy splitting between the "ground" states -E^o,e=±i(5') of the double- 
well Hamiltonian at the coupling parameter g = 0.005. In the first row, we indicate 
the quantity n which is the order of the expansion (j254p . In the second row, we list 
the partial sums (j254p . whereas in the third row, the Borel-Pade transformation is 
presented [cf. Eqs. (1220]) - (1222)) ]. 



n 


^ ^ 1^1=0(^0,1019 


rvnig) 









5.327 056 794 x 10" 


-14 








1 


5.169 464 698 x 10" 


-14 


5A78 046 952 x 


10" 


-14 


2 


5.166 551 926 x 10" 


-14 


5.166 468 039 x 


10" 


-14 


3 


5.166 379 090 x 10" 


-14 


5.166 367 389 x 


10" 


-14 


4 


5.166 364 871 x 10^^ 


-14 


5.166 363 195 x 


10" 


-14 


5 


5.166 363 450 x 10^^ 


-14 


5.166 363 250 x 


10" 


-14 


6 


5.166 363 286 x 10" 


-14 


5.166 363 262 x 


10" 


-14 


7 


5.166 363 265 x 10" 


-14 


5.166 363 261 x 


10" 


-14 


8 


5.166 363 262 x 10" 


-14 


5.166 363 261 x 


10" 


-14 


9 


5.166 363 261 x 10" 


-14 


5.166 363 261 x 


10" 


-14 


exact 


5.166 363 261 x 10" 


-14 


5.166 363 261 x 10" 


-14 



As in the case of the cubic anharmonic oscillator (j237p . these coefficients grow factorially and 
have the same sign. The perturbation series (j218p with these coefficients are, hence, divergent 
and not Borel-summable. As mentioned already in the discussion above, however, apart from 
the non-summability of the series (j218p another and even more critical argument against the 
application of the Rayleigh-Schrodinger perturbation theory to the double-well potential has been 
mentioned [178-181, 245-247]. Namely, while the perturbation series (j218p . if being resummed, 
would provide only one energy level for every quantum number A'^, two eigenstates are known to 
appear in the case of nonvanishing coupling g > 0. Such a splitting of the level of the unperturbed 
Hamiltonian into two levels (the instability of the eigenvalues from a mathematical point of view) 
has been explained by Landau and Lifshitz [237] within a semiclassical WKB approach and has 
been attributed to the effect of quantum tunneling through the potential barrier between the 
two wells. However, while the semiclassical approach has provided a qualitative and perhaps semi- 
quantitative understanding of the energy spectrum, it can hardly be extended for a full and accurate 
quantitative analysis of the double-well (as well as multi-well) potentials. 

In this direction, an alternative route for studying the eigenvalues of the double-well-like potentials 
has been proposed and described in detail in Refs. [178-181,245-247] and is based on the Feynman 
path integral approach. We restrict ourselves here to a rather short account of the basic formulas. 
In particular, it has been conjectured that the eigenvalues of the Hamiltonian (|247p can be found 
by solving the generalized Bohr-Sommerfeld quantization condition 

where e = ±1 is the parity, and the functions B(iw{E,g) and AdwiE,g) are determined by the 
perturbative expansion and the perturbative expansion about the instantons, respectively. The 
evaluation of these functions in terms of series in variables E and g has been described in detail 
elsewhere [179,180,247] for quite general classes of potentials. For the particular case of the double- 



AdAE,g) 



ei , 



(249) 
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well potential, for example, the function Bd^{E,g) has the following expansion: 



i?dw(i?, g) = E + g (?,E^ + +g^ {^^E^ + ^E^ + 0{g^) , (250) 

and defines the perturbation expansion (j218p which can be easily found by inverting the equation 
B^„{E, g) = N. The instanton contributions to the eigenvalues of the Fokker-Planck Hamiltonian 
are described by the function 

A^„{E,g) = ^ + g (UE^ + ^) +9' (227E' + ^e) + 0{g') . (251) 



By making use of the generalized Bohr-Sommerfeld quantization condition (j249p together with 
the expansions (|250p - (|25ip of the and B^w functions, we may finally derive the eigenvalues of 
the double-well Hamiltonian. That is, by expanding (|249p in powers of the two small parameters 
exp(— l/6g') and g, the eigenvalues of the Hamiltonian (j247p we may find a so-called generalized 
perturbation expansion (resurgent expansion) 

00 00 /i^\Nn / _1 /fir, \ " n-l 



/2\ / e~-^/^^\ f / 2\^'^ 

EnA9)=Y.^n9'+Y.[-) -"^^ E T.^^^r.ki9\ (252) 

1=0 n=i \ V y / k=o ^ \ y/ J 



which apart of the perturbation series also accounts for the instanton contributions. As seen from 
this expansion, apart from the principal quantum number A^, eigenvalues of the double-well po- 
tential are characterized also by the parity e = ±1. That is, each unperturbed level with the 
quantum number splits up into two states of opposite parity. Such a splitting is described by 
the second term in the right-hand side of the Eq. ()252p and is attributed to the instanton contri- 
butions. The order of the instanton contribution is defined by the index n, with n = 1 being the 
one-instanton contribution, etc. Despite of its name, the leading, one-instanton term involves a 
summation over all possible n-instanton configurations but neglects instanton interactions [179]. 
As seen from Eq. (1252p . the evaluation of the energies within such a (first-order) approximation, 
also referred as a "dilute instanton gas" approximation, requires the knowledge of the e^^nki coef- 
ficients. Since these coefficients are available for download [85], we recall here just the six leading 
ones for the ground state A^ = of the double-well potential: 

71 6299 2691107 

eo,ioo - 1 , eo,ioi - ~Y2 ' ~ ~"288" ' ~ ~ 10368 ' 

2125346615 509978166739 846134074443319 



60,104 .nv^^xi ' ^0,105 .0^71 n^c ' ^0,106 



497664 5971968 ' 429981696 

By inserting these coefficients in Eq. (I252p . we are able to perform a numerical check of the validity 
of the one-instanton expansion for the energy splitting 



A^o = ^0,-1 - ^0,1 ~ 2 eo,io; 9' (254) 



1=0 



of the ground-state level A^ = due to the quantum tunneling. As seen from the second column 
of the Tables [10] and [TT\ however, a straightforward, term-by-term summation of the series (j254p 
may well reproduce the energy splitting AEq for the case of the small coupling g = 0.005 but 
leads to the divergence of the result for strong coupling of g = 0.03. In order to overcome such a 
divergence we again employ the Borel-Pade resummation method as given by Eqs. (|220p - (|222p in 
which as input data we used coefficients eo,ioz- However, while the Borel-Pade method provides a 
good convergence of the series ()254p , the result of resummation shows significant discrepancy with 
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Table 11: Energy splitting between the "ground" states EQ^^=±i{g) of the double-well Hamiltonian 
at the coupling parameter g = 0.03. In the first row, we indicate the quantity n which is the order 
of the expansion (j254p . In the second row, we list the partial sums (j254p . whereas in the third row, 
the Borel~Pade transformation is presented [cf. Eqs. ()220p - (1222p ]. 



n 


^ ^ 2^1=0^0,1019 


TVnig) 









2.518 531 055 x 10"^ 








1 


2.071 491 792 x lO'^ 


2.178 799 750 x 


10" 


-2 


2 


2.021 916 083 x 10^^ 


2.011 031 238 X 


10" 


-2 


3 


2.004 265 988 x 10^^ 


1.989 410 405 x 


10" 


-2 


4 


1.995 553 827 x lO'^ 


1.983 547 787 x 


10" 


-2 


33 


- 1.487 012 669 x 10^ 


1.980 196 846 x 


10" 


-2 


34 


- 4.525 128 589 x 10^ 


1.980 196 848 x 


10" 


-2 


35 


- 1.418 424 114 X 10^ 


1.980 196 848 x 


10" 


-2 


36 


- 4.575 754 388 x 10^ 


1.980 196 848 x 


10" 


~2 


exact 


1.983 530 115 X 10"^ 


1.983 530 115 X 10" 


-2 



the "true" energy splitting as computed by the diagonalization of the Hamiltonian (j247p in the 
basis of the harmonic oscillator wavefunctions. Such a discrepancy indicates the importance of the 
higher-instanton terms which take into account the instanton interactions for the case of strong 
coupling. The evaluation of the higher-order corrections (n > 2) to the energies of the double-well 
potential is, however, a very difficult task since it requires a double resummation of the resurgent 
expansion, in powers of both g and exp{—l/6g). In this contribution, one may address the question 
whether it is possible to resum the divergent series that gives rise to the A^-w and B^w functions 
directly and look for solutions of the quantization condition (j249p without recourse to the resurgent 
expansion. As it was discussed recently [248], this appears to be the only feasible way to evaluate 
the multi-instanton expansion (in power of n) because the quantization condition incorporates all 
instanton orders. We briefly summarize the results of this investigation here, which has been carried 
out for a potential closely related to the symmetric double-well. 

We discuss the main ideas of such a "direct summation" approach for the case of the so-called 
Fokker-Planck potential: 

HFFig) = \v' + (1 - Vaqf + V9<i-\ 

= \p' + \q' -\ + ^q- V9q' + l'?'^ (255) 

which can be considered as a modification of the double-well potential (j247p . In contrast to the 
double-well case, however, the Fokker-Planck potential contains a linear symmetry-breaking term. 
Because of this term, the ground state of the potential (|255p is located in one of the wells and, hence, 
does not develop any degeneracy due to parity (see also the mathematically inspired discussion 
earlier in this Section of the review). Excited states, however, can still be degenerate for 17 — > 
in view of the (perturbatively broken) parity e = ±1 of the quantum eigenstates. Again, in order 
to analyze such a degeneracy we may start from the generalized Bohr-Sommerfeld quantization 
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condition which for the case of Fokker-Planck Hamiltonian reads as: 



Ti-BME,g))T{l-BFpiE,g)) \ gj 2ti 
where the series for the functions Bpp{E,g) and App{E,g) are: 



r,x 2BFpiE,g) -AFp{E,g) 

= 0, (256) 



BFp{E,g) = E + 3E^g + (^35E^ + ^E^ / + 0{g^) 



(257) 



AFp{E,g) = + ( 17i?" + 7 ) 5 + ( 227E^ + ^E]g' + 0{g-') . (258) 



and 



3g V 6 

As mentioned above, these expressions may be used directly in order to investigate the energy spec- 
trum of the Fokker-Planck potential. That is, if one defines the left-hand side of the quantization 
condition ()256p as a function of two variables E and g: 

The zeros of this function at some fixed g determine the energy spectrum of the ([255]): 

QiEN{g),g) =0, NeNo. (260) 

A numerical analysis of the function Q{E, g) can be used, therefore, in order to find the eigenvalues 
of the Fokker-Planck Hamiltonian. 

As seen from Eq. (|259p . an analysis of the function Q{E,g) can be traced back to the evaluation 
of the functions A-p-p{E,g) and B-p-p{E,g) which constitute the series in both variables E and g [cf. 
Eqs. (|257p - (|258p ]. In order to perform these summations, we may rewrite the functions A-p-p[E,g) 
and Bpp{E,g) as a (formal) power series in terms of some variable x, taken at x = 1: 



BME,g) =Y.bFkE,g)x%=i, (261) 

A:=0 



1 " 

App{E,g) = +Y,4kE,g)x%=i, (262) 
^ k=i 

where the coefficients a^p{E,g) and bpp{E,g) are defined uniquely from (j257p and (|258p : 
b^pl{E,g) = E, b^pl{E,g) = SE'^g, a^pl{E,g) = {iJE"^ + 5/6) g, etc. Of course, the power se- 
ries (|26ip and ()262p are mathematically equivalent to the expansions ()257p and ()258p but more 
convenient for further numerical computations as discussed below. 

Now we are ready to analyze the energy spectrum of the Fokker-Planck Hamiltonian [cf. Eq. (j260p ]. 
As mentioned above, to perform such an analysis for any particular g we have to (i) resum the series 
for the functions App{E,g) and Bpp{E,g) and (ii) to insert the resulting generalized Borel sums 
A-pp{E,g) and Bpp{E,g) into Eq. (|259p : these can be either the upper or the lower distributional 
Borel sums, but must be chosen consistently. We may interpret then the Q{E, g) as a function 
of E (at fixed g) and (iii) numerically determine the zeros of this function which correspond to 
the energy values -Eppf"* of the Fokker-Planck Hamiltonian [cf. Eq. ()260p ]. As an illustration of 
the calculations, we have computed the function Q{E,g) for the coupling g = 0.1 and the energy 
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Table 12: The function Q{E,g) [cf. Eq. ([259|) ] for the Fokker-Planck Hamiltonian 
calculated at the coupling g = 0.1 and the energy E = 5.199 138 696 x 10^'^. In ad- 
dition to this function (fourth column), the generalized Borel transforms Bpp{E,g) 
and App{E,g) of the series ()26ip and (|262p are displayed as a function of the order 
n (first column). Indeed, the exact values for the functions App{E,g) and B-pp{E,g) 
are not known. 



n App{E,g) Bpp{E,g) Q{E,g) 



1 


3, 


.421 


230 


+ 


0, 


,000 


000 i 


0, 


,005 


207 


+ 


0, 


,000 


000 i 


0.000 186 


+ 


0.000 175 i 


2 


3 


.418 


180 


+ 


0, 


.000 


000 i 


0, 


,005 


198 


+ 


0, 


,000 


000 i 


0.000 211 


+ 


0.000 175 i 


3 


3, 


.420 


640 


+ 


0, 


.000 


000 i 


0, 


,005 


372 


+ 


0, 


,000 


Olli 


0.000 031 


+ 


0.000 171 i 


4 


3 


.437 


970 


+ 


0, 


.015 


857 i 


0, 


,005 


366 


+ 


0, 


,000 


105 i 


-0.000 057 


+ 


0.000 006 i 


5 


3, 


.428 


560 


+ 


0, 


,019 


738 i 


0, 


,005 


366 


+ 


0, 


,000 


105 i 


-0.000 006 




0.000 026 i 


11 


3, 


.430 


080 


+ 


0, 


,017 


575 i 


0, 


,005 


355 


+ 


0, 


,000 


092 i 


-0.000 004 




0.000 002 i 


12 


3 


.429 


310 


+ 


0, 


,017 


691 i 


0, 


,005 


354 


+ 


0, 


,000 


090 i 


0.000 000 


+ 


0.000 000 i 


13 


3, 


.429 
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+ 


0, 


,017 


320 i 


0, 


,005 


354 


+ 


0, 


,000 


090 i 


0.000 000 


+ 


0.000 001 i 


14 


3 


.429 
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0, 


,017 


312i 


0, 


,005 
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0, 


,000 


090 i 


0.000 000 
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0.000 001 i 


15 


3, 
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+ 


0, 


.017 


321 i 


0, 


,005 


354 
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0, 


,000 


090 i 


0.000 000 


+ 


0.000 000 i 


16 


3, 


.429 


570 
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0, 


.017 


540 i 


0, 


,005 


354 


+ 


0, 


,000 


090 i 


0.000 000 


+ 


0.000 000 i 


17 


3 


.429 


470 


+ 


0, 


,017 


440 i 


0, 


,005 


354 


+ 


0, 


,000 


090 i 


0.000 000 


+ 


0.000 000 i 


18 


3 


.429 


470 


+ 


0, 


,017 


450 i 


0, 


,005 


354 


+ 


0, 


,000 


090 i 


0.000 000 


+ 


0.000 000 i 



exact 0.000 000 + 0.000 000 i 



E = 5.199 138 696 x 10~^. For this set of parameters, the series ()26ip - ()262p are divergent and we 
shall again use the Borel-Pade resummation method in order to obtain the functions A-pp{E,g) 
and Bpp{E,g). The resulting distributional Borel sums App{E,g) and Bpp{E,g) are presented 
in the second and third columns of Table [T2[ By inserting these generalized sums into Eq. (j259p 
the function Q{E,g) has been calculated and is displayed in the forth column. As seen from the 
Table, with the increasing the order n of the series (|26ip - (|262p . the function Q{E,g) approaches 
zero, indicating that the energy E = 5.199 138 696 x 10^^ is an eigenvalue of the Fokker-Planck 
Hamiltonian for a coupling g = 0.1. 

While the above calculations have been performed mainly for illustration purposes, the "direct 
summation" approach as given by Eqs. (|259p - (|262p has been recently employed for the detailed 
calculations of the ground as well as the excited state energies of the Fokker-Planck potential for 
a wide range of coupling parameters [248] . 



3.7 Renormalization Group and Critical Exponents 
3.7.1 Orientation 

Universal properties of critical phenomena have found a natural explanation thanks to the theory 
of renormalization group (RG, [249], [250]). This powerful conceptual tool provides a framework 
to understand the surprising behaviour that physical systems exhibit in the vicinity of second- 
order phase transitions (like universality, and the existence of a scaling behaviour parameterized 
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by critical exponents). Computations of critical quantities are intimately related to the techniques 
described in the present review: for instance, the so-called approach at fixed dimension d = 3 
combines the technical difficulties and sophistications of Feynman diagram calculations with the 
mathematical intricacies of the divergent character of the perturbative expansions; on the other 
hand a different method, based on so-called high-temperature series, leads to truncated convergent 
expansions whose analysis is connected to the numerical methods discussed in Sec. [2] of this review. 

However, before we analyze these two approaches in a more detailed manner, let us first recall some 
basic facts. Although this may perhaps come as a surprise, due to a number of technical reasons 
the extraction of high-precision estimates from the RG is a hard task, requiring a lot of effort and 
computational resources. As a matter of fact, even if much effort has been put in recent years 
towards a precise determination of critical quantities, only very recently have theoretical results 
reached an accuracy comparable to that of experiment; in addition, not always do the different 
methods necessarily agree among themselves, or with the experimental results. 

For example, the most precise theoretical estimate nowadays available (obtained in [251] from 
a sophisticated combination of Monte-Carlo and high-temperature techniques) predicts for the 
critical exponent a in systems described by an 0(2)-symmetric ((^^) -theory a value of a^"^ = 
—0.0151(3); this figures are incompatible with the currently accepted experimental result of a^p^ = 
—0.01285(38) (obtained on the Space Shuttle by the microgravity experiment MISTE [252,253]), 
and a convincing explanation for such a discrepancy still has to be found. On the other hand, the 
theoretical predictions produced by the class of the so-called field-theoretical perturbative methods 
appear to be in general agreement with the experiment, but usually at the price of a much larger 
uncertainty (giving for the quantity a^^^ mentioned above estimates like a^r^ = —0.011(4) of 
Ref. [254], where the theoretical error has been estimated conservatively, or the less conservative 
a^=2 = -0.0129(6) of Ref. [255]). 

In any case much more precise experiments which could help to solve these open questions are 
planned for the future, even if the recent changes in NASA's priorities have unfortunately lead to 
the cancellation of an entire family of devices which were almost ready to fly (see, e.g., Ref. [256] 
for their detailed description). 

In general, many computational tools have been developed to formulate numerical predictions 
for universal physical quantities like critical exponents: among them we find lattice-based methods 
(such as Monte-Carlo methods [257-261] and high-temperature series [262,263]), perturbative field- 
theoretical methods (like the e-expansion [264-271], the approach at fixed dimension d = 3 (see 
Refs. [272-274]) or the minimal renormalization without e-expansion [275]), and non-perturbative 
field-theoretical methods (based on the so-called "exact" RGs, see, e.g., Refs. [276-279]). Almost 
all of these techniques are based on the construction of more and more precise approximations for 
the desired quantities, which naturally lead to the evaluation of a truncated series. 

In Refs. [264-266], the so-called e-expansion led to the first convincing a priori prediction for 
critical exponents, thus validating the RG approach. A significant breakthrough toward the precise 
computation of RG quantities has been achieved soon after by another perturbative field-theoretical 
method, the approach at fixed dimension d = 3, originally suggested in [280,281]. In a seminal series 
of articles [272-274], the calculation of the RG /3-function for the 0(A^)-symmetric (i/)^) -theory 
in d = 3 advanced from the 2-loop level (see [280,282]) to the level of 6 loops. The treatment of 
Refs. [272-274] (consisting in the numerical computation of ~ 1000 integrals in dimensions up to 
6 with 8-10 significant digits) opened the possibility of formulating quite precise predictions for a 
large number of physical quantities, which as of today remain among the best ones available. Quite 
recently, other methods (notably Monte-Carlo estimates, eventually supplemented by the output 
of some ad hoc high-temperature series) have reached a comparable or even better accuracy (see 
for example [251]). 
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It would be very difficult (and out of the scope of the current review) even only to mention all 
the existing methods to compute RG-related quantities. For example, the complete description 
of the techniques employed to obtain critical quantities by means of just one single method (the 
5-loop evaluation of the e-expansion) is of the size of a book [283]. Even more so, a whole series of 
books [284] has been devoted over the years to keep track of all progress in the field. 

As mentioned before, in the present work we will try to give an overview of two methods in 
particular, (i) high-temperature series (see Sec. 13.7.2^ and (ii) the approach at fixed dimension 
d = 3 (see Sec. 13.7.3]) . The choice of focusing on these two approaches to critical systems is justified 
by the following considerations. 

Both methods deliver their results in form of series; however, reflecting the very different nature 
of the two computational schemes the estimates derived from HT series usually possess a finite 
radius of convergence, while the approach at finite dimension produces descriptions of the critical 
quantities in terms of divergent series; in addition, the number of accessible terms of the series is 
very different in the two cases (typically ~ 25 for the HT expansions, and only 6-7 terms for the 
scheme at fixed dimension). 

From a practical point of view, however, many analogies between the two methods can also be 
found: in both cases the computation of the coefficients of the series requires the evaluation of 
a factorially-growing number of different contributions in higher order; in both cases adequate 
prescriptions (an acceleration scheme for the first method, and a resummation for the second) 
have to be supplemented if one wishes to be able to extract meaningful physical results from the 
obtained series. Quite surprisingly, the final precision provided by the many terms of HT series is 
usually comparable to that obtained from the few terms of the scheme in fixed dimension. 

Indeed, the two methods show an intriguing mixture of similarities and differences, both being 
good examples of situations where a sophisticated numerical analysis of some series is helpful or 
mandatory in obtaining refined physical results; for this reason we think that a detailed inspection 
and comparison of such methods should be considered particularly interesting in the framework of 
the present review. 

Finally, in subsection 13.7.41 we present a general outlook of the perspectives to extend the results 
obtained from both the HT series expansions and the techniques at fixed dimension. 

3.7.2 High— Temperature Series 

The general idea of high-temperature expansions may be understood as follows. Let us assume that 
we are dealing with a statistical-mechanics system which undergoes some (second-order) phase 
transition and is described by some suitable partition function Z, which depends on the inverse 
temperature /? = \/{kBT). An approximate description of the system valid only at temperatures 
T ^ ^critical niay then be derived if we expand the partition function in terms of /? ~ l/T. (Of 
course, a complementary "low temperature" -expansion valid when T <^ To-iticai may also be derived; 
the two expansions are related, and a way of exploiting such connection is described in [285]). 

An interesting and non-trivial property of such an expansion is the possibility of writing it in purely 
diagrammatic terms; even if the number of diagrams which must be generated and evaluated grows 
very quickly with the perturbation order, only sophisticated graph-generation and combinatorial 
techniques are required to complete the evaluation of the series. This feature should be contrasted 
with the method at fixed dimension explained in the following subsection, where not only the 
number of diagrams, but also the computational complexity of the integral associated to each 
diagram grows with the order in perturbation theory. 

Many evaluation schemes for HT series have been presented in the literature (for an introduction 
to this subject the reader is referred to [286]); a very successful specialization of the technique is 



78 



given by the so-called linked-cluster expansion (LCE) as described in [287]. We usually find the 
LCE applied to a class of physical systems described by the following partition function, which is 
a generalization of the Ising model: 



e/5E(.,) _ (263) 



Here, (pi is a scalar field, while / is a non-negative distribution function which decreases faster than 
e""^ as 4> —>■ oo and is normalized by the relation f dcj) f{4>) = 1; the sum runs over all pairs of 
nearest neighbors on a lattice which is usually taken to be a generalization of quadratic lattices to 
arbitrary dimensions. 

An interesting feature of this way of setting up the problem is the possibility of encapsulating the 
physical content of many different systems into the choice of the distribution function / which 
enters Eq. (I263p : for example, the choice 




/W= >^ 5|* + 2»y-j^^j (264) 

reproduces the usual spin-S Ising model, while many other models which are well-known in statis- 
tical mechanics (like Blume-Capel, Klauder, double Gaussian and range models) can be recovered 
by different choices. In addition, the computation of the LCE can be set up in such a way that its 
final results are directly expressed in terms of the symbolic expressions for the cumulant moments 
of /. These quantities ^2n are defined by the relation 

E(£;T^'" = l-/d'/'/('A)e''^ (265) 

Thus, once the result is obtained as a function of the ^2n) it can be specialized to any of the possible 
models just by replacing the symbolic moments with the actual values they assume when a specific 
choice of / has been made. 

The physically interesting quantities to be computed are as usual the connected (/-point functions 
at zero magnetic field defined by 



9" \nZ[h] 



dhi^ . . . dhig 



(266) 

h=0 



with their moments 



M^"^ = Yl {xi-X2fGq{x^,...,Xq). (267) 



Xq = 

X2,...,Xq 

r(l) 



X2,...,Xq 



The sophisticated graph-theoretical methods which allow the high-temperature series to be actually 
generated are described, e.g., in Ref. [288]. We just mention in passing that, once the form of the 
lattice and the expansion variables have been chosen, the core of the method is the possibility of 
rewriting the truncated expansion of the sum on the right-hand side of Eq. (1263^ as the sum of 
many contributions, each one associated to a particular graph; the new summation runs over a 
set of graphs satisfying some topological constraint, and it is possible to formulate graphical rules 
to describe the coefficient attached to each graph of the expansion. In addition, many different 
expansions may be generated, due to the fact that the contributions to the original expansion can in 
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turn be rearranged in terms of new contributions associated to a different set of graphs. In the end, 
the method consists in the (recursive) generation of all the graphs satisfying the constraints dictated 
by the chosen expansion, together with the numbers specifying their associated contributions. 

Here, we would like to discuss in some more detail how the obtained data can be interpreted. A 
typical HT series reads (see [288]) 

X2 = l + 8?; + 56?;^ + 392t;3 + 2648v^ + 17864w^ + 118760v^ + 789032z;^ 

+ 5201048?;^ +34268104 + 2246798647;^° + 1472595144 1;^^ +9619740648?;^^ 
+ 62823141192 z;^^ + 409297617672 + 2665987056200 + 17333875251192 -y^^ 
+ 112680746646856 v^"^ + 731466943653464 z;^^ + 4747546469665832 1;^^ 
+ 30779106675700312 7;20 _^ 199518218638233896 1;^^ + 1292141318087690824 1;^^ 

+ 8367300424426139624 1;^^ + 54141252229349325768?;^'' 

+ 350288350314921653160 v"^^ + 0(^2^) , (268) 

with V = tanh (3. We notice that the generated series can be quite long (the typical order be- 
ing ~ 25), and that they have been computed for many different physical systems and lattice 
arrangements; an extensive list of such computations can be found for example in [289]. 

However, despite their length HT series usually show poor convergence, and sophisticated methods 
must be employed to analyze their content in terms of RG quantities. A detailed account of many 
of the techniques which may be used (not only in relation to the case of HT series, but also for 
generic similar statistical- mechanical problems) can be found in [30]. The typical model problem, 
as derived from the predictions of the RG theory or statistical mechanics and considered in this 
reference, is how to fit the obtained HT series S{z) = Yl^=Q^nZ^ to a function of the following 
form: 



S{z) ~ A 1 



l + Ai[l--] \a2(i--\ ' + ... 



(269) 



for z ^ Zc from below, with < Ai < A2 < .... The constants A, Zc and A are called the 
critical amplitude, the critical point and the critical exponent, respectively. The terms of the form 
(1 — z/ Zc)~^\ which may be present or not, and eventually be replaced by logarithmic singularities 
like |ln(l — z/zc)\^\ are called corrections to scaling, and if one or more of the (or the 5j) are 
non-integral we say we are in the presence of confluent singularities. It is perhaps interesting to 
remark that for some special systems even more complicated forms for the critical singularities are 
actually predicted by RG theory, and observed experimentally. 

As shown in [30], the efficient extraction of the values for A, Zc and A (plus eventually the Ai 
and the Aj) from series like Eq. (j268p is far from being straightforward. The most frequently used 
approaches fall in two main broad classes: 

(i) Methods consisting in fitting/extrapolating the sequence of the ratios of the terms of the 
series, that is the new series Rn ^ Sn/Sn-i- The main reason for employing the ratio series 
lies in the fact that it has better convergence properties. In addition, under reasonable as- 
sumptions about the functions A and Ai appearing in Eq. (|269p . it is possible to rewrite the 
ratio terms as (see [30]) 



5„, 1 



A — 1 c d e 

+ + ^;T+K + ;^TT2A + ;^TT3A+ (270) 



Sn~l ^critical 

+ 0(n^2) ^ o(n^2-A) ^ o(n^2-2A) ^ o(n"2-3A^ ^ _ 

which is a form particularly suitable for a fit versus 1/n. (Actually, the last formula is valid in 
the case of one confluent singularity, but it may be generalized to more complicated instances. 
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The easiness of the fit can be readily checked for the simpler case of no confluent singularity, 
which is recovered by setting A = 1.) 

(ii) Methods derived from Pade approximants. Notably, a generalization of Fade approximants 
called differential approximants — which has been proposed in [290] and subsequently devel- 
oped in Refs. [291,292] and other works — has been derived just for the specific purpose of 
extracting the values of critical quantities, and stands nowadays as the technique of choice 
for the analysis of series like Eq. (j268p . 

Here we will content ourselves with a few observations more about the results obtained by using the 
methods included in class (i) above, to make a connection to Sec. [2j Interestingly, many if not all of 
the "classical" acceleration schemes presented in Sec. I2.2.71 (Richardson extrapolation, e-algorithm, 
0-algorithm, Levin transforms, and so on) have been applied in the literature to the study of ratio 
series of the type of Eq. (1270p . The careful and detailed analysis presented in [30] — carried out 
comparing many different methods both for some artificial test series and for some series derived 
from real problems in statistical mechanics — shows how some of the examined acceleration schemes 
(notably the ^-algorithm) behave more favourably than others, consistently with the observation 
that for the case of ratio methods a logarithmic convergence should be expected (see below); even 
more consistently, methods mainly suitable to accelerating linear convergence, like the e-algorithm, 
score poorly in such a situation. 

One may now ask how it is possible to show that the convergence of the ratio sequence Rn should 
be logarithmic. The answer lies in the fact that the terms of a ratio series can be seen as originating 
from a partial summation of a series Pn defined as 

D _ E> D Sn-l SnSn-2 — (Sn-l) /o'7i ^ 

= -Kn - — T; T; — n 5 • ' J- J 

Jn-l 'Jn-2 'Jn~lJn~2 

Assuming the general form (j270p for the ratios, it can be checked that 

constant 



Pr. 



^critical 



A-1 c{l + A) ^ 



n(n-l) n(n-l) 



i+A 



11^ 



(272) 



so that Pn/Pn-i ~^ 1; ^ud wc havc logarithmic convergence indeed. 

This example is instructive since it shows that there are real-life situations where the validity of 
Eq. (f27l) can be directly tested. 



3.7.3 Approach at Fixed Dimension d = 3 

In QED, the fine-structure constant as the fundamental coupling parameters is small, and the 
computation of more terms in the perturbative expansions is usually effective to get more precise 
estimates of physical quantities even if the series we start with would not in principle be convergent. 
However, in the case of the determination of critical quantities by means of many field-theoretical 
methods, the addition of terms of the expansion derived from higher loops is simply not enough 
even when it is feasible: the coupling constant is so large in the relevant region where the phase 
transition occurs that the additional terms do not by themselves provide any increase in the 
precision obtained for the results. 

This consideration explains why in the latter case the loop calculations must be complemented by 
a resummation prescription in order to determine critical exponents. In the case of the method 
at fixed dimension d = 3, this observation is even more legitimate, since in such a framework the 
Borel summability of scalar 0{N) theories has been rigorously proven [293]. In this light, it is not 
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any more surprising that from an historical point of view, the computation of critical exponents by 
field-theoretical methods has represented a testbed for the application of many new resummation 
schemes; in particular, many of the devices implying Borel summation have been introduced in 
physics or invented just with the purpose of producing better estimates for critical exponents. 

Typical results obtained for the /3- function and the anomalous dimensions r]{g) and 772(5) hi the 
0(2)-symmetric three-dimensional (</)^) ^-theory read [274] 

-u + u"^ - 0.4029629630 + 0.314916942 

-0.31792848 + 0.3911025 - 0.552448 u'^ + 0{u^) , (273a) 

0.0118518519 - 0.0009873601 + 0.0018368107 
-0.00058633 + 0.0012514 + 0{u^) , (273b) 

2 4 

— u + —u'^- 0.0495134445 + 0.0407881056 
5 50 

-0.04376196 + 0.0555573 +0{u'^). (273c) 

Analogous series have been obtained for other values of N. Following RG theory, a different uni- 
versality class (that is, a different set of physical systems) corresponds to each value of this pa- 
rameter. For example, the case = 2 is found to describe superfluid He'^ at the A-point and 
Xy-ferromagnets, the case = 1 labels Ising-like ferromagnets and liquid-to- vapour transitions, 
the case A = 5 has been proposed as descriptive of some high-Tc superconductors [294], and the 
limit A — > models polymers and self-avoiding walks [295]. 

These series describe the coefficients of the Callan-Symanzik equation 

A + f3(g) ^ _ _ ft) ^(g) - 6^2(5)) r(-''') = Ar(-''') « (274) 

for the so-called three-dimensional 0(A^)-vector model, which is defined by the bare Euclidean 
action (the index B labels the bare quantities) 

5b(03) ^ / Q(V03)2 + + f (^1)^) . (275) 

Here, <j)^ is a vector of A scalar fields, F'-'^''') are the one-particle irreducible Green functions 
(correlators) for the renormalized theory, with a external vertices and b insertions of 0^ operators. 
The quantities g and u are related by a change of normalization. We recall that the right-hand 
side of Eq. (j274p depends in general on the renormalization scheme adopted [see again [274] for 
the precise definition of a scheme leading to a right-hand side as in (I274p ]. 

We would like to emphasize once again that the numerical coefficients in Eqs. (|273ap — (|273cp are 
deceptively simple quantities. In fact, to obtain just each last term of these series, on the order of 
1000 Feynman amplitudes had to be integrated numerically, with a precision of 8 — 10 digits each. 
Such amplitudes were originally known in terms of integral representations in 18 dimensions, and 
only a cleverly arranged set of computational simplifications allowed for a reformulation [274] in 
terms of equivalent amplitudes of lower dimensionality (< 6), thus opening the way to a successful 
high-precision numerical evaluation. 

The goal of calculating critical quantities in this computational scheme is not straightforward, even 
if the coefficients in Eq. (j273p are assumed to be computable numerically. Until the end of this 
subsection, we examine in detail all the necessary steps to be performed. 

(i) The diagrammatic expansion. First of all, we have to generate all the diagrammatic expansion 
for the Green functions r(^'°), F^^'") and F^^'^), together with that of dV^"^'^^ / dp^ . The knowledge 
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Incremental number 


Total number 


Loop number 





1 


2 


3 4 5 6 7 


8 


6 7 8 


r{2.o) 








1 


2 6 19 75 317 


1622 


103 420 2042 


r{4.o) 


1 


1 


2 


8 27 129 660 3986 


26540 


828 4814 31354 



Table 13: Number of diagrams (one-particle irreducible and without tadpoles) con- 
tributing to the two-point and four-point function up to 8 loops in perturbative 
order. 



of these correlators is indispensable to have enough information to be able to extract RG functions 
P, r] and 7/2 from Eq. (j274|) . respecting at the same time the renormalization prescriptions set by 
the renormalization scheme employed. 

A measure of the complexity of the generation stage can be found in Table [T3l It should be noticed 
that the number of diagrams for QF^^'*') /dp"^ is excluded from the Table, the corresponding raw 
number of diagrams being empirically intermediate between that of F^^'*^) and that of F^^'*^^. On 
the other hand, since our renormalization scheme is such that the correlators must be evaluated at 
zero external momenta, the diagrams of F(2.i) can be considered a subset of those of F^^'*^). 

For each amplitude, the corresponding combinatorial/symmetry factor, and the corresponding 
0{N) polynomial must be computed too. 

(ii) Finding the best parameterizations for the needed amplitudes. At a given loop level L, the natural 
representation of each amplitude entails a 3-L-dimensional momentum-space integral. Although 
other representations are in principle possible, it was realized in [274] that, starting from the 
standard momentum representation of diagrams, one can systematically find and replace in the 
amplitudes one-loop functions by analytic expressions which are known in d = 3 dimensions. This 
procedure lowers the number of needed final integrations. 

As an example, we will show in some detail how the parameterization of the diagram in Fig. [5][^a) 
can be carried out. 

The diagram in Fig. [5ja) has four loops. After redrawing it as in Fig. E^b) to highlight its planar 
structure, and assigning loop momenta as in Fig. [5]^c), we can make use of standard Feynman rules 
for scalar theories in 3 dimensions, arriving at the formula 

Ao^ j d^^i d3^2 d^^3 d^^4 Pr^(^i) Pr(4) Pr(£4) 

X Pr(£i + ^2) Pr(£2 + t-^) Pr(£2 + ^4) Pr(4 - t-i) , (276) 

which expresses the amplitude at zero external momenta (as specified by the renormalization 
scheme adopted in [281] and [274]) in terms of an integral representation of dimensionality 12. The 
propagator Pr(€) here is the scalar Euclidean kernel 

£^ -|- 

The renormalized mass m is set to one by the renormalization schemes used in [281] and [274]. 
Needless to say, starting from Eq. (j276p . it would be extremely difficult, if not impossible, to give 
a precise numerical estimate for the amplitude corresponding to the diagram in Fig. [5]^a). 

Now, we can formulate a much more efficient strategy to evaluate the diagram in Fig. O^a), if we 
suppose for example to know, with its complete dependence on external momenta ei and 62, the 



83 




(a) 






(h) 



(i) 



Figure 5: Parameterization of a four-loop diagram in terms of "triangles" [diagrams 
(a), (b) and (c)], and various steps in the parameterization of the diagram in Fig.[5|^a) 
in terms of effective triangle vertices [diagrams (d)-(i)]. Further explanations are 
given in the text. 



analytic formula for the one-loop 3-dimensional integral corresponding to the "triangle" diagram 
of Fig. El^d), that is, 

Tri(ei,e2) = j d^£ Pr(£) Pr(£-Fei) Pr(£ + 61-^62) . (278) 

In fact, in this case we can identify in the right-hand side of Eq. (j276p suitable subintegrals being 
of the form of the right-hand side of Eq. (j278p . replacing them with the analytic formula for the 
function Tri; this operation can even be performed diagrammatically, by systematically substituting 
all occurrences of Fig. [5][[d) with the "effective vertex" of Fig. EJ^e), where the dashed lines attached 
to triangle vertices are used to automatically implement momentum conservation. 

This method leads, after the simple graphical manipulations of Figs, ^g), (h) and the choice of 
momenta of Fig. EJ^i), to the new integral representation 

^oc / d^£i d^£2 Tri(£i,-£i) Tri(£i,£2) Pr(£i+£2) Pr(^2) • (279) 
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After the explicit parameterization of momentum space d^ii d^£2, this representation turns out 
to have a dimensionahty of "only" 3. In Ref. [273] a value of 0.1120344277(1) is quoted for this 
diagram (we notice that a suitable normalization factor therein described must be inserted into 
Eq. ()279p to obtain this numerical result.) 

As the four-loop diagram in Fig. E^a) constitutes a relatively simple example, we would like to 
add that much more complex effective vertices than the triangle alone must be used to make 
this technique effective for the large sets of diagrams which are needed in the six- and seven-loop 
orders. In addition, the way of replacing effective vertices is not unique, possibly leading to many 
different integral representations, each one with a different final dimensionality. For example, as 
Fig. [5|^f) already suggests, many other triangles than the two highlighted ones could be identified 
at different positions in the diagram (however, it should also be noticed that in this example the 
maximal number of triangles which can be replaced at the same time is always two, no matter 
which choice of triangles is considered) . The identification of the most favorable parameterization 
defines a very complicated optimization problem, which, if connected to the very large number of 
diagrams to be evaluated, calls in turn for an automated way of finding out the most favorable 
parameterizations. A strategy to solve this problem in a more general and detailed way, even 
considering other possible analytic simplifications, is described in [296]. 

It should be noticed that in this approach the renormalization procedure can be built into the step 
of amplitude parameterization. The reason is that the O (N) vector model is superrenormalizable 
in 3 dimensions (see [296] for further details). 

(iii) Numerical evaluation of the needed amplitudes. Each integral representation obtained dur- 
ing the last stage must be numerically evaluated at zero external momenta, following the 
renormalization-scheme prescriptions originally formulated in Ref. [281]. When combined with the 
symmetry factors and the O (N) factors, as described above, this operation leads to the perturba- 
tive knowledge (up to the desired loop order) of RG functions /3, rj and r]2 appearing in Eq. (j274p . 

As previously emphasized, the difficulty implied by this stage should not be underestimated. For 
example, no published algorithm for multidimensional numerical integration is able to guarantee 
that the needed number of significant digits (7-8) will be obtained after a reasonable number of 
function evaluations for all needed 9-dimensional integral representations at the level of 7 perturba- 
tive loops. However, an innovative algorithm described in [297], which has been recently designed 
with this problem in mind, should finally allow to obtain an asymptotic behaviour of the integration 
error nice enough to reach the high precision required for this multi-dimensional computation. 

(iv) Resummation of the obtained series. In order to complete the analysis of the critical exponents, 
we recall that the non-trivial value of the coupling u such that the RG /3-function in Eq. (I273ap 
goes to zero must be identified numerically (this corresponds to the so-called Wilson-Fisher fixed 
point of the O (N) vector model). One can then extract, for example, the value for the critical 
exponent a from the following formulas 

/3{u*) = 0, (280a) 

« = , . . (280b) 

2 — r][u*) + r]2{u*) 

In this physical situation, the critical value of the coupling is non-perturbatively large (we have 
g* ~ 1.4 here), and therefore a suitable resummation of the series for the scaling functions is 
required to get meaningful results for the desired physical quantities. In addition, and as mentioned 
before, the mathematical property of Borel summability has been rigorously established in [293] 
for the ((^^) theory in the case d = 3, so applying the procedure of Borel resummation to this 
physical problem is definitely legitimate. 

From an historical perspective, it is noteworthy that the first precise numerical evaluation of critical 
quantities has been derived in [274] by means of a Borel-Pade summation as described in Sec. 13.2.31 



85 



However, a spectacular theoretical breakthrough was achieved in the seminal paper [176,298] when 
some additional information, derived from the large-order behavior analysis of perturbation theory 
for scalar quantum field theories carried out in [299], has been conveyed to the Borel resummation 
of these series (a detailed review of the subject with additional results can also be found in [300]). 
Such analysis states that the coefficients of the perturbative series for any critical quantity S{g) 
behave at large perturbative orders n as 

5„ = c(-a)Vn! (^1 - O (^^^^ , (281) 

with a, b and c being computable numbers which depend on the critical function to be evaluated. 
In particular, a turns out to be independent of the series under consideration and equal to 

9 

a = 0.147774232(1^ 



iV + 8 



for the case d = 3. It should be noticed how this discovery legitimates a posteriori the use of the 
Borel-Pade summation in [274], which had been proposed in that context as a mere guess. 

If we introduce a generalized Borel-Leroy transform b9^'^ \u) of S{g), defined as 



n=0 

(which differs from the formulas presented in Sec. 13.2.11 only for the presence of an arbitrary 
parameter b'), Eq. ()28ip shows that the closest singularity of ' (n) (called instanton singularity) 
is located at the point — 1/a (we refer the reader to [177] for a more detailed discussion about both 
this fact and the used terminology). 

If the additional hypothesis is formulated that all the finite-distance singularities of the Borel 
transform of the series come from subleading instantonic contributions (which are located further 
away on the negative axis of the complex plane), and the Borel transform is thus supposed to be 
analytic in the whole complex plane except for the cut going from —1/a to — oo (this statement 
is also known as the maximal analyticity hypothesis), then a conformal mapping specially tailored 
to the analytic structure of the series can be used, leading to a remarkable enhancement in the 
accuracy of the obtained resummed quantities. It should be noticed that, contrary to the property 
of Borel summability, the property of maximal analyticity has not been rigorously proven for the 
case of d = 3; however, arguments for its plausibility are given in [177]. 

We then arrive again at a method based on the conformal mapping of the Borel plane, which we 
have already discussed in detail in Sec. I3.2.4t by introducing the transformation 



^/TTau- 1 

z [u) = , (283) 

y/l + au + 1 

we map the Borel plane, cut at the instanton singularity u = —1/a, onto a circle in the z plane 
in such a way to enforce maximal analyticity, and thus to enhance the rate of convergence. With 
the help of Eq. (j283p . we can now rewrite the Borel transform ()282p upon expansion of the u in 
powers of z: 

oo 

B^^'''\u{z)) = Y^ 4''!^^ (284) 

n=0 
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where the coefficients Sn ' can be uniquely determined as function of Sn- FinaUy, the fuh physical 
solution S{g) can be reconstructed by evaluating the Laplace-Borel integral: 



oo 




n=0 



Interestingly, it is possible to connect the so far arbitrary parameter b' to b, and use it as a 
variational check to validate the method; in addition, introducing more complicated prescriptions 
the knowledge of c can be incorporated in the resummation algorithm as well (the reader interested 
in this point may wish to consult [176]). In this way, one can properly take advantage of all the 
available information about the analytical properties of the critical quantities. Finally, additional 
refinements of the method described here are possible, like considering an additional homographic 
transformation to displace possible complex singularities in the (/-plane [254]. 

However, we would like to mention that many other resummation methods have been proposed and 
implemented in this context. For example, an alternative approach which does not make use neither 
of Borel summation nor of Fade approximant is that of order-dependent mapping given in [210] 
and [211]; to the best of our knowledge, this method too has been invented for the particular case 
of the resummation of critical exponents. 

Anyway, let alone the problem of finding the correct value for the resummed series, we would like 
to emphasize how an even greater care should be spent in giving an accurate evaluation of the final 
uncertainty of the resummation. In fact, the error in critical quantities basically comes from two 
sources, the summation error at fixed g* , AS, and the error induced by the error in g* , Ag*: 



or the series 1/S, can give precious hints about the spread of the numerical results, thus ultimately 
leading to a reasonable indication about the summation errors. A sophisticated analysis of the 
error bounds, as well as a careful review of the values for critical quantities obtained by Borel 
resummation from quantum field theoretical techniques can be found in [254] (see also the earlier 
work [298]). 

3.7.4 Concluding Remarks on Renormalization Group Related Quantities 

Since the determination of critical exponents is a particularly prominent application of numerical 
convergence "accelerators" and summation techniques for divergent series, we would like to include 
here a few remarks regarding possible foreseeable progress in this field. Indeed, there is some hope 
that the computations presented in the two last subsections can be partly extended in a not too 
distant future. 

As discussed in [301], in the case of the HT series the main computational problem is apparently a 
graph-theoretical one. In particular, the algorithm of generation of the set of graph required for the 




(286) 



In this respect, the computation of some auxiliary quantities like the shifted series 




(287) 
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LCE is very inefficient, being plagued by the so-called problem of isomorphic copies: every known 
algorithm for graph enumeration generates some subset of the needed graphs more than once, so 
that the effective number of graphs to be scanned during the generation stage can be orders of 
magnitude larger than the dimension of the desired final set. In addition, storage problems are 
quite common during the computation of the series (see Refs. [259-261]). As a matter of fact, the 
reached level of 20-25 terms of HT series seems to be close to the saturation of presently available 
computational resources. 

Similar considerations can be carried out in the case of the approach at fixed dimension: even if 
promising steps have been done in [296] towards the full extension to 7 perturbative loops of the 
original 6-loop computation of [274], the method does not appear to be easily extensible to a higher 
perturbative level. The most relevant problem with such an approach lies in the fact that there is 
no a priori proof of its validity: no theorem exists stating that at a given loop order, for a given set 
of available effective vertices, some minimal reduction in the number of residual integrations will be 
obtained uniformly for all graphs needed for the evaluation of the desired physical quantities; the 
effectiveness of the method can in fact be demonstrated only a posteriori, by an explicit inspection 
carried out diagram by diagram. Indeed, such a proof has been recently given in [296] for the level 
of 7 perturbative loops; however, it is unlikely that the method can be extended much further, 
since each loop more raises the upper bound for the dimensionality of the residual integration 
which must be carried out for each Feynman diagram, thus quickly making things computationally 
too expensive even if the considered set of effective vertices is quite large. 

Other computational schemes for RG quantities which are less suffering from this explosive growth 
in computational complexity with the order of the approximation do exist — notably Monte-Carlo 
methods. As mentioned before, it is possible with these methods to reach a precision comparable 
to that of the experimental results so far available; however, also the scaling of these algorithms 
is actually very bad, requiring as well a formidable consumption of computational resources (for 
example, the sophisticated theoretical estimates quoted in [251] required about 20 years of total 
CPU time). 

In any case, while some progress can certainly be achieved with all three of the above mentioned 
methods, it appears unlikely that any of the theoretical techniques known so far will be able to reach 
a precision comparable to that foreseen for the next generation of experiments (which should be 
at least one order of magnitude better than the accuracy reached so far) . Perhaps some theoretical 
breakthrough could change this picture. 

4 Conclusions 

First, let us give a synopsis of the convergence acceleration methods discussed in Sec. [2] of this 
review. We can summarize a few important aspects as follows: 

(i) As has been emphasized in [31] and here, our theoretical knowledge of the properties of many 
of the algorithms is incomplete, although an important classification can be given in terms 
of alternating and nonalternating, as well as linearly and logarithmically convergent series. 
We can say that a theory of their properties, if it is to be considered as satisfactory from 
a practical point of view, should provide a theoretical estimate of the numerical stability 
properties in addition to the convergence acceleration properties of the transforms. 

(ii) None of the convergence accelerators discussed here is linear in the sense that the nth-order 
transform of the sum of two input series is equal to the sum of the nth-order transforms of 
the two input series. Linearity is thus not provided. However, a weaker requirement called 
translation under translation can in general be imposed [see, e.g., Eq. (3.1-4) of Ref. [31]]. It is 
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generally acknowledged nowadays that the nonlinear transformations discussed in this review 
potentially lead to much better numerical results that typical linear convergence acceleration 
schemes. 

Logarithmically convergent sequences are rather difficult to accelerate; in particular, they 
constitute a set with the property of remanence (see Ref. [25]), thus satisfying the hypothesis 
of a "non- universality" theorem [79]. This implies that no universal sequence transformation 
exists which is able to accelerate all the sequences belonging to the set; in other words, a 
given method will work only for some subset of the logarithmic sequences, and will not work 
for all the others. However, a few potentially very powerful algorithms can be indicated (see, 
e.g., Sec. [2X2] and Table [I]). 

It would be a hopeless task to attempt an exhaustive overview of all available convergence 
acceleration methods in this review. As an example of an algorithm which is not discussed 
here, we only mention the Overholt process and the ACCES algorithms which are described 
in detail in Ref. [25], or the Euler transformation and the Euler-Maclaurin formula which 
are discussed in many standard reference books on mathematical methods for the physics 
community (e.g., Ref. [76, 161]). We have only attempted to provide a general overview of 
some available methods, to introduce a rather universally applicable construction principle 
via the iJ-algorithm (see Sec. I2.2.7p . and to illustrate these general concepts by considering 
some examples. The interested reader will find further information in the more specialized 
literature indicated in Sec. [2j 

If one is simply interested in numerical results, there may be a strong motivation to skip the 
theoretical phase of the algorithm-finding process sketched in Sec. [2] to a large extent, and 
to investigate the performance of convergence accelerators experimentally. It is undeniable 
that this approach if far from satisfactory from a theoretical point of view, but it can lead, 
nevertheless, to numerically spectacular results. In the case of no available analytic infor- 
mation, experiments are even unavoidable. In this context, one might note that according 
to Refs. [106,302], the nonlinear sequence transformation listed in Table [H notably the d 
and 6 transformations discussed in Sec. 12.2.31 below, are among the most powerful and most 
versatile sequence transformations for alternating series that are currently known. 

The notion of a certain "complementarity" appears to be rather universally applicable to 
the field of convergence acceleration algorithms, in various contexts. E.g., the e algorithm 
is known to be an efficient accelerator for linearly convergent and even factorially divergent 
input series, but fails for logarithmically convergent input data. The p algorithm is known 
to be an efficient accelerator for logarithmically convergent input series, but fails for linearly 
convergent and factorially divergent input data. Both the e and p algorithms are numerically 
rather stable and rather universally applicable to series in their respective domains, but are 
in many cases outperformed by special sequence transformations which profit from explicit 
remainder estimates. On the other hand, the sequence transformation with remainder fails 
completely if the assumptions about the structure of the remainder term are not met by a 
given input series. 

The notion of "complementarity" also affects the "iteratability" of the transformations: Rel- 
atively simple sequence transformations like the process, which make only very crude 
assumptions about the structure of the remainder term, can typically be iterated, leading to 
improvements in higher orders of the iteration process. In sharp contrast, the elimination of 
a complex remainder structure as implied by a sequence transformation with a remainder 
estimate intuitively leads to the conjecture that a further iteration of the process cannot 
increase the apparent convergence any further. This is indeed confirmed in numerical studies 
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performed by the authors, and illustrates the general complementarity aspect of the conver- 
gence accelerators once more: the power of a specialized algorithm has to be traded for a 
disadvantage in the sense that further iterations typically cannot increase the numerical gain 
of a single iteration. 

The perhaps most important application of the convergence acceleration techniques is the calcu- 
lation of the hydrogen Green function, described in Sec. 12.4.21 Namely, while the so-called Za- 
expansion ( "perturbative bound-state QED") in general converges better for low nuclear charges 
Z than for high Z, it is necessary, at the current level of accuracy, to even include nonperturbative 
effects at low Z . This surprising observation is especially important for the one-loop self-energy 
correction, which is the dominant radiative correction in hydrogenlike ions. To give an example, the 
difference between all known terms of the Za-expansion up to the order a{Zafmc^ (Ref. [303]) 
and the nonperturbative, higher-order remainder for the IS state [35] amounts to 28 kHz. This is 
about 700 times larger than the current best measurement in atomic hydrogen, which has an ac- 
curacy of about 40 Hz [304,305]. Thus, while corrections of smaller absolute magnitude can still be 
described perturbatively at low Z (e.g., two-photon effects), the dominant one-photon corrections 
have to be calculated nonperturbatively to match current experimental precision, and from this 
point of view, atomic hydrogen (Z = 1) constitutes a "high-Z system" whose description would 
be incomplete without a calculation of the nonperturbative effects. An essential ingredient in that 
direction is given by the methods discussed in Sec. I2.4.2[ 

Regarding divergent series and Borel summation treated in Sec. [3l one important conclusion is 
that even the best summation method cannot reconstruct a complete physical answer from per- 
turbation theory alone, if there are instanton-related contributions, nonperturbatively suppressed, 
which must be added to the perturbative expansion. This has been seen in one-dimensional model 
problems, notably in the case of double-well like potentials treated in Sec. 13.61 In the calcula- 
tion [248], an attempt has been made to use resurgent expansions at moderate and large couplings, 
for a few nontrivial cases of one-dimensional potentials which, nevertheless, admit a treatment via 
generalized quantization conditions. The conclusion has been that while the resurgent expansion 
can be useful to explore nonperturbative effects at small coupling, it cannot be used for accu- 
rate quantitative methods in the asymptotic strong-coupling region. Roughly speaking, the triple 
resummation of the resurgent series in n, and / in Eq. (j252p is too hard for any available sum- 
mation algorithms to handle even if one enters directly into the quantization condition ()249p and 
directly resums the A^^ and -B^w function given in Eqs. (j250p and ()25ip . Nevertheless, the accurate 
understanding of nonperturbative effects in quantum potentials with instanton contributions via 
resurgent expansions at small and moderate coupling could be seen as a significant consequence 
of the understanding gained by the path integral formalism. Let us also mention that important 
connections of the quantum oscillators to field theoretical problems exist, as explained in detail in 
Ref. [177]. 

We should perhaps mention that our treatment of the even and odd anharmonic oscillators in 
Sees. 13.41 and 13.51 is entirely based on the concept of distributional Borel summation and thus far 
from exhaustive. One aspect left out in this review concerns renormalized strong coupling expan- 
sions [306,307] and related methods (see also [102,108,111,112,307-310]) which make possible the 
determination of very accurate numerical energy eigenvalues using only few perturbative coeffi- 
cients, and which cover both small and large couplings. Expansions of this kind, however, fail for 
double-well like potentials and also cannot reproduce the width of resonances in the case of an odd 
perturbation. 

To conclude this review, let us indicate some open problems and directions for further research, 
(i) There is a regrettable absence of a general mathematical proof for the convergence of some 
of the nonlinear sequence transformations described in Sec. [2j This might lead to a direction of 
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research for mathematics. Only for some special model problems could rigorous convergence proofs 
be obtained (see [311-313] or Sections 13 and 14 of [31]). (ii) New sequence transformations with 
improved remainder estimates could be a significant direction of research for both physics and 
mathematics. In this context, one might mention hyperasymptotics (see Ref. [76] and a couple of 
preparatory calculations, e.g., Refs. [314,315]. (iii) For nonalternating divergent series, the next 
term omitted from the partial sum of the asymptotic series is typically not a good simple remainder 
estimate. One may indicate the unsolved problem of finding a practically useful remainder estimates 
for nonalternating divergent series, in particular those which give rise to an imaginary part in 
the sense of a distributional Borel summation, (iv) Let us also indicate a possibly interesting 
further direction in the field-theoretical direction which is based on the following observation: In 
general, a perturbation potential or a field interaction determines the large-order behaviour of the 
generated perturbation series. However, the interaction potential also determines the large-coupling 
expansion of the related quantities. This means that there should be a connection between the large- 
coupling expansion and the large-order behaviour of the perturbation series. This notion and related 
ideas about connections of weak- and strong-coupling regimes have quite recently been explored in 
Refs. [66-69,255,316,317], but further work in this direction is definitely needed, (v) Finally, there 
is a possible connection of the "prediction" or extrapolation of perturbative coefficients as described 
here in Sec. 12.2.51 and Brodsky-Lepage-Mackenzie (ELM) scale setting [318] which is based on 
the idea that, roughly speaking, seemingly wildly divergent perturbative expansions in field theory 
can sometimes be formulated as not so wildly divergent series if the expansion variable (coupling 
constant) is taken at a momentum scale appropriate to the problem. This approach is illustrated 
in Ref. [318] using the two-loop correction to the muon anomalous magnetic as an example, and 
a related approach has recently been generalized to higher orders in Ref. [319]. However, it seems 
that there is still room for improvement in the renormalization-group-inspired approaches to the 
extrapolation of perturbative expansions to higher-order loop corrections in field theory. 
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A Basic Terminology for the Spectra of Hamilton Operators 

A.l Resolvent, Spectrum, Green's Function 

In this subsection we shortly recall, for the convenience of the reader, some of the basic mathe- 
matical definitions concerning spectral theory of Schrodinger operators, including resonances. A 
standard reference for these topics is [44]. 

Definition 16 Given a closed operator H in a Hilhert space 7i with domain D{H), its resolvent 
set p{H) is the set of all X G <C such that the operator H — XI (where I is the identity operator) is 
invertible and its range R{H — XI) is the whole ofH. 
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By the closed graph theorem (see e.g. [320], Theorem III. 5. 20) the inverse operator R{\,H) = 
{H — XI)~^, which by the above definition exists and has domain Ti, is bounded, or, equivalently, 
continuous, namely there is C(A) > such that, for all u £ TC: 

\\R{X,H)u\\ <C{X)\\u\\. (288) 

p{H) is called resolvent set because the linear inhomogeneous equation (H — XI)u = / in the 
unknown u £ Tl can always be solved: it has indeed the unique solution u = R{X,H)f for any 
datum / G W if A e p{H). 

Definition 17 The complement of p{H) with respect to <C, i.e. the set of all X G C such that either 
R{X,H) does not exist (in the sense that the operator H — XI is not invertible), or, if it does, its 
domain is not 7i, is the spectrum of H, denoted cr{H). 

The following examples can be given: 

1. Let H he a n X n matrix acting on 7i = C". Its spectrum cr{H) is clearly the set of the 
eigenvalues of H, because H — XI is not invertible if and only if Hu = Xu for some u ^ 0, i.e. 
det(-ff — XI) = 0. On the other hand, li H — XI is invertible, its range is always 7i because 
7i has finite dimension. 

2. Likewise, in general any eigenvalue E oi H belongs to cr{II), because for any corresponding 
eigenvector ip one has {H—EI)'ip = 0. Therefore the spectrum always contains the eigenvalues. 

3. If the Hilbert space is infinite dimensional, however, the spectrum may have a totally different 
nature. Let for instance H = X he the position operator in quantum mechanics, namely the 
maximal multiplication operator by x in L^(1R,). Its action and domain are: 

(Xn)(x) = xu{x) , D{X) = {u£ L'^{]R) : xu{x) £ L^{U)} = 

u{x) : / |u(x)p dx < +oo; / \xu{x)\'^ dx < +oo [ . (289) 
Ju Jn J 

Then a{X) = R. We obviously have, indeed: 



({R(X,X)u){x) = ^. (290) 
X — X 



Now, for any given A G R, pick u{x) G L^(1R,) continuous and nonzero in an open interval 



around A. Then ["^^L dx does not exist and hence R{X, A) is not defined on the whole of 



n. 

The most important example in the present context is represented by the Schrodinger operator 
H = T + V acting in ^^(IR,'^). Here T is the kinetic energy operator, and V the (real-valued) 
potential. In units with h = 1, these read 

(rn)(x) = -An(x) , u £ D{T) , {Vu){x) = V{x)u{x) , u £ D{V) . (291) 

Here, D{T) and D{V) are the maximal domains of T and V, namely 

D{T) = {u£ L2(r") : Au £ ^^(R")}, D{V) = {u £ ^^(R") : Vu £ l2(R")} . (292) 

We take as the domain D{H) the maximal one, 

D{H) = {u£ L2(R") : -Au + Vu£ ^^(R")} . (293) 
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It is easy to see that T and V are self-adjoint operators, namely T = T* , V = V* , where A* 
denotes the adjoint of the operator A. Under suitable conditions on V (see Ref. [44], Chap. 10), 
also H is self-adjoint (with the proof being often far from trivial). Concerning self-adjointness, we 
recall that: 

1. If A is any self-adjoint operator in Ti., then (t{A) C R. 

2. The symmetry property {Au\v) = {u\Av), V(n, f) G D{A) is not sufficient to guarantee 
self-adjointness because it implies only A*u = Au for u £ D{A) but it does not exclude 
D{A) 7^ D{A*). The difference is quite substantial because while the spectrum of a self- 
adjoint operator is real, the spectrum of a symmetric but not self-adjoint operator is the 
whole of C 

By its very nature, any Schrodinger operator is a differential operator acting on a Hilbert space 
of functions. This makes it possible to characterize its resolvent and spectrum in terms of the 
solutions of the time-independent Schrodinger equation {—A + V)'ip = E'ip. 



Definition 18 For any A G C, let Q{x] A) he a distributional solution of the inhomogeneous equa- 
tion 

-Ag{x;X) + [V{x)-X\g{x;X)=6{x), x G R" , (294) 

with 6{x) being the Dirac distribution. Under very general conditions on V , it can be proven that 
Q{x;X) is unique, and it is called fundamental solution, or Green function, of the Schrodinger 
operator H = T + V . 



The Green function yields an explicit representation of the resolvent [H — A/]~^. Consider indeed 
the integral operator in L^(1R,"') of kernel Q{x — y;X), namely the operator acting on functions 
u £ in the following way: 

(G(A) u) {x)= f Q{x- y- A) u{y) dy . (295) 
We have, for any u £ 1? (the interchange between differentiation and integration can be justified): 
(i/-A/)(G(A)n)(x) = [-A + y(x)-A] / Q[x - y-\)u{y)dy = (296) 



[-A + y(x)-A]6;(x-y;A)n(y)dy= / b{x - y) u{y) dy = u{x) . (297) 

R" JR" 

This shows that G(A) is the inverse of \H — \1\, and therefore coincides with the resolvent operator. 
Hence the above formula yields the explicit action of the resolvent -R(A, -ff): 

RiX,H)u{x)= [ g{x-y-\)u{y)dy. (298) 

JR" 

As an example, consider the free particle energy operator in one dimension, 

H = T = DiT) = {u:u£ L^,u' £ L^,u" £ LH . (299) 

dx^ 

The differential equation for the Green function is therefore 

-g"{x]\)-\g{x-\) = 5{x). (300) 
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To solve this equation, we switch to Fourier space. Let g{p; A) be the Fourier transform of g{x; A): 

g{p;X) = ^ [ g{x;X)e-'P''dx. (301) 
v27r Ju 



Then (integrating twice by parts) we have 

^ / (x; A) e-'P" dx = p^g{p; A) . (302) 
v2vr Ju 



27r Jk, 

Since J^6{x)e~^P^ dx = 1, the differential equation becomes 

= ^ e,,;,)=_L_^ (303) 

If we now invert the Fourier transform (using for instance complex integration) we readily obtain 

g(x;A) =e^l^l. (304) 

Hence, we have 

R{T, X)u =[T- XI]^^ u{x) = [ e^l^-^l u{y) dy . (305) 

From this representation, we immediately see that cr{T) = [0, +oo) as it should be. If indeed 
A E C, A ^ ]R,+ , then we have for all u £ L'^, with A = |A| exp(i(/)): 



\\R{T,X)u\\^< / / |e2^(^-^)| 1^(^)12 dydx< / \u{y)\'^ dy / ^inW2))\x-y\ 

a/|A| sm(0/2) Jr ^|A|sm(0/2) 

whence the boundedness of R{X,T) because < C(A) < +oo as long as < cp < 2-k, i.e. A ^ 
[0,+oo). In three dimensions with x = (xi, X2, X3), in exactly the same way we obtain for the 
Green function of T = — A: 

^\\x\\ 



g(x;A) = ^^, \\x\\ = ^^xl + xl + xl. (307) 
As above, we can check that (t{T) = [0, +cxd). 

In the above examples, the spectrum is continuous. The question emerges how to characterize the 
spectrum in terms of the solutions of the Schrodinger equation. We recall that by definition, E is 
an eigenvalue of the Schrodinger operator H (and hence a quantum bound state) if there is a non- 
zero function il}{x,E) E L^(1R"') (and hence normalizable: in other words, \iIj{x^E)\^ dx = 1) 
such that i^V' = Eip. Then, 'ijj{x,E) is called an eigenfunction corresponding to E. The number 
m of linearly independent solutions of Hip = Eip is the multiplicity of the eigenvalue E. 
By the normalization condition, all eigenfunctions corresponding to any bound state vanish as 
||x|| — > +00. Hence E is an eigenvalue of the Schrodinger operator if and only if the Schrodinger 
equation Hip = Eip admits at least one normalizable solution. It is well known that eigenfunctions 
corresponding to different eigenvalues of a symmetric operator are orthogonal: 

{ip{x,E)\iP{x,E')) = I iP{x,E)lp{x,E')dx = 0, E ^ E' . (308) 

Linearly independent eigenfunctions corresponding to the same eigenvalue can always be chosen to 
be orthogonal by the Gram-Schmidt orthogonalization procedure when 7i is separable, i.e. when 7i 
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contains a countable basis. It follows that the eigenfunctions corresponding to the bound states form 
an orthonormal system. If the spectrum consists entirely of eigenvalues, this system is complete, 
namely every vector in TC can be expanded in Fourier series of the eigenvectors. 

Under very general conditions on V, it can actually be proven that any eigenfunction is continuous 
and vanishes at least exponentially (pointwise), i.e. there is K[E) > such that 

— > oo . (309) 

The points belonging to the continuous spectrum of H can be characterized in terms of the solutions 
of the stationary Schrodinger equation Hip = having a totally different behaviour at infinity: 

A point A G R belongs to the continuous spectrum of H if and only if there is a solution f{x, A) ^ 
of the Schrodinger equation H f = Xf which fulfills the pointwise bound \f{x, A)| < K(X) for some 
K{X) > and all x G R". 

Again, it is useful to consider an example. Consider once again the kinetic energy T in one dimen- 
sion. Setting A = fc^, /c G R, the Schrodinger equation Tij) = Eip becomes 

- -r^i^ix, k) = k^Mx, k) (310) 
dx^ 

with the well known plane wave solutions 

i;{x,k) = e^'''^ , (311) 

which are not in because \ip{x,k)\ = 1 if and only if A; G R; however ijj{x,k) is obviously 
bounded if and only if fc G R, and hence A > 0. The functions ip{x, k) are often called continuum 
eigenfunctions or sometimes generalized eigenfunctions. Since j^e^^^ dx = 2ir6{k), where 6{k) is 
the Dirac distribution, they fulfill the continuum orthogonality relations 

[ i){x,k)^{x,k')dx = 5{k-k'). (312) 

Sometimes, this relation is written under the form {'ip{x,k)\tp{x,k')) = 6{k — k') even if the scalar 
product in the right-hand side does not make sense because ipix, k) ^ L^. 

Solutions of the Schrodinger equation which are both non-normalizable and unbounded have no 
relation with the spectrum. 

We recall that an orthonormal system (ON) in a (separable) Hilbert space is a (finite or count- 
able) set of vectors ij^n-, with n G No, such that (V'nlV'm) = ^m,n- Given any vector n G its 
Fourier coefficients with respect to the given ON system are by definition the scalar products 

On = (ti|^n) • (313) 
The ON system is complete when any vector in 7i can be expanded in Fourier series, namely when 

oo 

U = ^ On'tpn , (314) 
n=0 

where the convergence takes place in the norm of TC. 

If the spectrum of the Schrodinger operator H is discrete, namely it contains only isolated eigen- 
values of finite multiplicity, then the set of all its eigenfunctions can be selected in such a way that 
it forms a complete ON system. The same property holds in general for any self-adjoint operator 
with discrete spectrum. 
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Let us consider again some example. Let H be the quantum harmonic oscillator, H = —-^ + x'^ 
defined on its maximal domain. Then its spectrum is discrete, because it can be proven that it 
contains only the simple eigenvalues En = 2n+l with n G Nq. The corresponding eigenfunctions are 
the Hermite functions i^n{x) = , ^ ^ e~^' Hn{x), where Hn{x) is the nth Hermite polynomial. 

y 2" n! y/if 

The system ipn is orthonormal and complete: 

u{x) = 'S^ an-ipnix), an = —= ^ / u{x) e'"^ Hn{x) dx . (315) 

The convergence takes place in the norm (or quadratic mean), namely the above equality sign 
means 

N 

lim / \u{x) — } ttn ipn{x)\'^ dx = , (316) 

N^OO ha 

^ n.=0 

and neither implies pointwise convergence nor uniform convergence. 

If the spectrum contains a continuous component, the ON set formed by the eigenfunctions is not 
complete even when it is countable. For example, the eigenfunctions of the Hydrogen atom are not 
complete in L^(IR^), because the spectrum of the corresponding Schrodinger operator H = — A — ^ 
has a continuous component along [0, +oo) in addition to the bound states En = 

When the spectrum contains a continuous component the completeness statement requires the 
consideration of all eigenfunctions, namely also the continuum eigenfunctions, which are the only 
ones when the spectrum is purely continuous. 

Assume for the sake of simplicity that the continuous component stretches over [0,+oo), so that 
setting k = \J —E we can label the corresponding eigenfunctions by (;/)(x, /c) = e^^^' /^/2tt for A; S R. 
Then the continuum Fourier coefficient of u{x) G LP' is by definition the integral transform 



c{k) = / u{x)(t){x,k)dx . (317) 

If, as above, we denote by tpn{x) the eigenfunctions corresponding to the discrete eigenvalues En, 
the completeness theorem for the eigenfunctions is expressed by the formula 

oo „ 

u{x) = an ipn{x) + / c{k) (f>{x , k) dk . (318) 

n=0 

Two examples illustrate the above statements, (i) Consider once again the kinetic energy T in 
one dimension. There are no eigenvalues, and as we know the spectrum is purely continuous along 
[0, +(X)). Then, the completeness theorem provides nothing else than the Fourier transform. We 
have indeed 

c{k) = [ u{x) e-''^^ dx = u{k) , (319) 

v27r Jr 

where 'u{k) is the Fourier transform of u. The above completeness formula is here just the Fourier 
inversion formula 

u{x) = [ c{k) (j){x, k) dk = [ u{k) e'^^ dk , (320) 

V27r V2vr Jr 



and the orthogonality relation ^^il:[x,k)Tp{x,k')dx = 6{k — k') immediately yields the Fourier- 
Plancherel formula: 

/ \u{x)\'^dx= [ \u{k)\'^dk. (321) 

(ii) Consider the radial equation for the Hydrogen atom, H = —A — -p-. Then Eq. (I318P holds 
for any function /(r) G L^(1R_|_) with ipn{x) the Laguerre functions of argument ^ and (j){r,k) a 
hypergeometric function as defined, e.g. in [237]. 
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A. 2 Analytic Continuation of the Resolvent and Resonances 



Quantum mechanical resonances are usually defined by complex, second-sheet poles of the scat- 
tering amplitude. This means what follows: one first proves that the scattering amplitude is an 
analytic function of the energy in the whole complex plane cut along the positive real axis. Then, 
if the potential is regular enough, one proves also that the scattering amplitude admits a meromor- 
phic continuation across the cut. The complex poles of the continuation are the resonances. The 
real part of the pole is the position of the resonance, and the imaginary part the width. It can be 
shown, under fairly general conditions, that it suffices to prove these properties for the resolvent 
of the Schrodinger operator. Namely, one proves by the Lippmann-Schwinger equation (see, e.g., 
Ref. [321]) that the existence of complex poles of the analytic continuation of the resolvent (shortly, 
second sheet poles of the resolvent) implies that the scattering amplitude admits the very same 
second sheet poles. 

Let us examine more closely the second sheet poles of the resolvent. First recall that, given any 
closed operator H, the resolvent R{X, H) is an analytic operator-valued function of A € p{H) 
(shortly, the resolvent is analytic for A G p{H)). This means that for any ip G TC the scalar product 

R4X) = {iP\R{X,H)ij) (322) 

is an analytic function of A G p{H). If H is self-adjoint, R^,{X) is in general analytic for A ^ R. If 
(t{H) is discrete, R^{X) is meromorphic with simple poles along the real axis at the eigenvalues of 
H. Let us illustrate these situations by means of two examples from quantum mechanics. 

(i) The position operator X. We have already seen that R(^X, is just the multiplication by . 
Hence: ^ 

R{X,X)i^{x) = =^ {^\R{X,X)ij) = [ ^-^i^^^dx, A^R. (323) 

X- X X - X 

The above formula for i?^(A) = {ip\R{X, X)tp) in general defines two distinct analytic functions, 
one in the upper half-plane ImA > denoted by R^{X) and the other one, denoted R^{X), in the 
lower half-plane Im A < 0. They are not the analytic continuation of each other, unless '0 vanishes 
in some open interval. Moreover, the possibility of continuing i?^(A) to the lower half-plane or 

R^{X) to the upper half-plane depends on the regularity properties of Consider for instance the 
case il^{x) = exp(— x^/2), and let us prove that i?^(A), a priori analytic for ImA > 0, admits an 
analytic continuation to the whole half-plane ImA < 0. In the same way, R^{X), a priori analytic 
for ImA < 0, admits an analytic continuation to the whole half-plane ImA > 0. We have 

i?+(A) = / - — -dx, ImA> 0. (324) 

Since the integrand is analytic in x for Imx < and vanishes faster than exponentially when 
|Rea;| — > -|-oo, we can shift the integration down by —a, namely we can perform the change of 
variable x — > x — ia, a > 0. Then 

KW = iz—r-,'i':- (325) 

The integrand in the right-hand side is analytic in A for ImA > —a, and the integral converges 
absolutely and uniformly for the same values of A. Hence it represents the analytic continuation of 
R^{X) to the half-plane ImA > —a. Since a is arbitrary, iZ^(A) thus admits a continuation to the 
whole of C A symmetric argument holds for i?^(A). 

(ii) For the harmonic oscillator, let us verify directly that if H is the harmonic oscillator, the 
function R^{X) is, for any ^ G L^, meromorphic with simple poles, which may occur only at the 
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points at A,„ = 2n + 1 with n S No- Recall that if ^„ is a normalized eigenvector of H corresponding 
to the eigenvalue = 2n + 1, then 

R{\,H)lPn = [H-XI]'^iJn = ^^^Pn, (326) 

An — A 

whence 

Let now tp ^ L'^ he arbitrary. By Eqs. (j314|) and (|315p . we can write 

oo 

■0 = ^ a„ ^pn, an = (V'lV'n) , (328) 

n=0 

whence 

oo oo oo 

R{X, H)ij = R{X, H) «n, ^n = Y.an RiX, H)^n = Y. ^\ _ ^ ipn , (329a) 

n=0 n=0 n=0 

SO that 

oo 

R^X) = {^IAR{X,H)^) = ^ (a^V'^l—^^Vn) 

m,n=0 

oo _ oo I |2 

^ 2n + l- X^^^^ ' ^2n + l-X ^ ' 

m,n=0 n=0 

The last equality follows by the orthonormality relation (V'mlV'n) = Sm,n- Hence IR^(A) is a mero- 
morphic function with simple poles which may occur only at the points A = 2n + 1 with n £ No- 
Notice that the pole at A = 2A; + 1 is present if and only if the Fourier coefficient = (V'lV'fc) is 
non-zero. 

We can now give the definition of resonance associated with a given Schrodinger operator H = T+V 
having continuous spectrum. For the sake of the simplicity we will admit that the continuous 
spectrum stretches either from zero to +00 (which is always the case when y — > as ||x|| — > +00) 
or from —00 to +00 (as in the Stark effect). 

Definition 19 A complex number E = Ei + iE2, E2 < 0, is a resonance of the Schrodinger 
operator H as above if there exists a dense set Q inTC such that 

1. For any ^ G Q, the function R^{X) = {'ip\R{X, H)ip) , a priori holomorphic for ImA > 0, 
admits a meromorphic continuation to ImA < across the cut generated by the continuous 
spectrum; 

2. E is a simple pole of any continuation R^{X). 

We remark that the assumption £^2 < is standard and corresponds formally to a resonant state 
being a decaying state in the future. Consider indeed the time-dependent Schrodinger equation 
Hil) = idtip. Assume that there exists 0(2;) not in which however is a pointwise solution of 
the time-independent Schrodinger equation Hcj) = E(f) with E as above. Then (^(x, t) = </>(x)e~'^* 
solves Hip = idtip and |(/>(x, t)\ = \(j){x)\e^^^ — > exponentially fast for any x as t — > +00 if E2 < 0. 
This is also the reason why the resonance width is the inverse life-time of the resonant state. 
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In proving the existence of resonances for a given Schrodinger operator the above definition requires 
the solution of two mathematical problems. The first is represented by the actual construction of 
the meromorphic continuation of the matrix elements i?,/,(A) across the cut for any ip belonging 
to a suitable dense set in TC. The second is the proof that the continuation actually has poles. 
The most powerful method to construct the continuation is the dilation analyticity of Balslev and 
Combes [191] (and its variants). It requires the analyticity of the potential V, and works also for 
the Stark effect [219]. Let us shortly review how it works. 

Given -0 G L^(1R,^), and 9 £ the corresponding dilated vector is 



^,(x) = e^''/V(e'x) . 
Notice that ipg and tp have the same norm in because 



IIV'.(^)II^ 



e X 



' dx 



|V'(x)|2dx 



(330) 



(331) 



The dilated position is clearly e^x, and the dilated momentum e ^Vx- Hence the dilated Schro- 
dinger equation is 

H{e) = -e-2^A + V{e^x) . (332) 
If we add a linear potential F xi, with F > 0, we get the Schrodinger equation of the Stark effect 

Hf = -A + V{x) + Fxi, (333) 

whose dilated form is 



Hf{9) = -e^20A + Vie^x) + e'^ F xi . 
Let now 9 be complex, < Im^ < 7r/4. The following statements can be proven: 



(334) 



1. The complex- dilated operator H[9) has the same real eigenvalues as H. For < argA < 
2arg0, the spectrum of H{9) consists at most of eigenvalues, which do not depend on 9. 

2. The complex dilated operator Hf{9) has discrete spectrum, i.e. it consists at most of isolated 
eigenvalues with finite multiplicity, with negative imaginary part. 

3. Let ifj £ LP' he a dilation analytic vector. This means that the function {ipQ, cp) is an analytic 
function of 9 for |Im ^| < 7r/4 for any (p £ L?' . The set of the dilation analytic vectors is 
dense in L^. Then one has: 



1 



1 



Hf-XI 



HF{e)-xi 



ipe 



|Im 9\<^, 



O>Im0>-f . 



(335a) 
(335b) 



The above formulae yield the explicit analytic continuations of the matrix elements R^{H) = 
{'ip\ [H — \I]~^'ip) and R^{Hp) = {tp\ [Hp — XI]~^tp), respectively. Moreover, we see that the poles 
of the continuation must coincide with the complex eigenvalues of H{9) and Hp{9), respectively, 
which are independent of 9. This last characterization of the resonances, namely as eigenvalues of 
the non self-adjoint operator H{9) uniquely associated with the given operator H by the complex 
dilation is particularly useful. 

The above result settles the question of the definition of resonances and of the construction of the 
analytic continuation of the resolvents (when the potentials are analytic). Concerning the actual 
existence of resonances, their characterization in terms of complex eigenvalues of the complex 
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dilated operators H{9) and Hf{9) has made it possible to obtain a rigorous existence proof in the 
physically most relevant cases (atomic systems, Stark effect [321,322]). 

Finally, we remark that the resonance E = Ei-\-\E2 is an eigenvalue of the complex dilated operator 
H{9). Therefore it corresponds to an eigenvector ^/^(x, E, 6), which depends on ^. If — > 0, ip{x, E, 9) 
becomes a solution of the stationary Schrodinger equation, H ■ip{x, E) = Eip{x, E), which is neither 
nor oscillatory because E is complex. In one dimension typically it exhibits exponential increase 
in one direction and exponential decrease in the opposite one. 

B Distributional Borel Summability 

B.l Mathematical Foundations of Distributional Borel Summability 

In Sec. 13. 2. n we have already seen that the Borel method cannot be directly applied to the divergent 
expansion Yl'n'=o d"^ because its Borel transform has a pole on the real positive axis and thus 
prevents the convergence of the Borel-Laplace integral [see the discussion following Eq. (jl78p ]. 
The same consideration applies to any divergent expansion F{g) ~ X^^o 9^ where the an ~ 
r(A;n + l) display a nonalternating sign pattern. Without loss of generality, we may assume a„ > 
for all n. The Borel transform B^p\g), defined according to Eq. (I176p . has a positive radius of 
convergence. However, since the coefficients of its expansion are all positive, we can state that 
by the Vivanti-Pringsheim theorem [323-325], the function B^p\g) has a singularity at the point 
where the convergence circle meets the positive real axis. 

As we have seen in Sec. 13.2.11 the series Yl'^^=o 9^ should be associated to the function M{g) = 
dte~^/{l — gt) for g > 0, where the integral can either be interpreted as a Cauchy principal 
value integral or as an integral encircling the pole at t = 1/g in one of the possible directions. 
Intuitively, one might conjecture that at least the result of the principal-value integration, given 
this particular integration prescription, should be uniquely associated to the divergent input series, 
whereas the sign of the imaginary part is ambiguous. The main problem is to determine the 
uniqueness conditions of the association of the divergent input series to the value of the principal- 
value integral, which in turn may be referred to as a "distributional" Borel sum, for reasons to be 
discussed below. The problem of clarifying the association is solved by the method of distributional 
Borel summability [189]. Essentially, this concept extends the Nevanlinna criterion to boundary 
values of analytic functions, which are in general distributions, whence the name. The formal 
definition is as follows: 

Definition 20 We say that the formal power series f{g) ~ Yl^=o^"-9^ Borel summable in the 
distributional sense to f{g) for < g < R, with R > 0, if the following conditions are satisfied: 

(i) In accordance with Eq. Ijil73^ , we set 



Then B{t), which is a priori holomorphic in some circle \t\ < A, admits a holomorphic 
continuation to the intersection of some neighbourhood of R+ = € IR : t > 0} with 
C+ = {t € C : Imt > 0}. 

(ii) The boundary value distributions i?(tibiO) = lim^^o -^(tibi e) existMt € 1R,+ and the following 
representation holds, for g belonging to the Nevanlinna disk Cr = {g : Iieg~^ > R~^}, 



oo 




(336) 




(337) 
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where PP(5(t)) = I {B{t + iO) + B{t + iO)). 



As for the ordinary Borel sum, the representation (|337p is unique among all real functions admitting 
the prescribed formal power series expansion and fulfilling suitable analyticity requirements and 
remainder estimates. In the ordinary Borel sum case these are the Nevanlinna conditions of Theorem 
m Their generalization to the distributional case is as given in the following Theorem [189]. 



Theorem 21 Let f{g) he hounded and analytic in the Nevanlinna disk Cr 
and let f{g) = i {4>{g) + 'Pig)), with (j)(g) analytic in Cr so that 



m - E 



On 5 



n=0 



< Aa{e)^ N\\g 



\N 



ViV= 1,2, 



(338) 



uniformly in Cr^^ = {g & Cr : aigg > —tt/2 + e}, Ve > 0. Then the series Yl'^=oi^n/n-\)u^ is 
convergent for small \u\, and its sum admits an analytic continuation B{u) = Bi{u) + B2{u). Here, 
Bi{u) is analytic in = {u : |Imn| < c?^^}, and B2{u) is analytic in Cj = {n : (Imu > 0, Keu > 
—d~^) or \u\ < d^^} for some d > 0. B{u) satisfies 



A' ft 
\B{t + \r]Q)\ < — exp - 



(339) 



uniformly for t > 0, for any rjQ such that < rjo < d ^. Moreover, setting PP(i?(t)) = ^ {B{t + 
iO) + B{t + iO)), the distributional Borel sum f{g) admits the integral representation: 



f{g) = - rvV{B{t)) exp 
9 Jo 



--] dt. 

9. 



9^Cr, 



(340) 



i.e. f{g) is the distributional Borel sum of Y1'^=q o-n 9^ forO < g < R in the sense of Definition [gQl 
Conversely, if B{u) = ^^=Q{cLn/n\)u"' is convergent for \u\ < d~^ and admits the decomposition 
B{u) = Bi{u) + B2{u) with the above quoted properties, then the function defined by ^340^ is 
real-analytic in Cr and (j){g) = B{t + iO) exp {—t/g) dt is analytic and satisfies i338\) in Cr. 



Proof See [189]. 

It is convenient to include the following remark. The function (^{g) = g~^ /q°° -B(i+iO) exjp{—t/g) dt 
is called the "upper sum" and (j){g) = g~^ B{t + iO) exp{—t/g) dt the "lower sum" of the series 
X^^o ^" 9^- follows that, for g > 0, f{g) = Re (j){g). On the other hand, with this method, we can 
single out a unique function with zero asymptotic power series expansion, that is the "discontinuity" 

d{9) = - r \B{t + iO) -B{t + iO)] exp ( --) dt = 4>{9) -W) ■ (341) 
9 Jo ^ J V 5/ 

Thus, d{g) = 2ilm (j){g), for g > 0. 

Important example cases, where the concept of distributional Borel summability finds applications, 
are being discussed in Sees. [331 (odd anharmonic oscillator) and below in Sec. IB.2I (Stark effect). 

Exactly as for the ordinary Borel summability [we recall Eq. p77|) and Theorem [B], the distribu- 
tional case can also be extended to perform the summation of power series diverging as fast as n!*^, 
A; G N. The corresponding generalization of Definition [20] is: 
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Definition 22 Let k be a rational number. We say that the formal power series F{g) ~ X^n^o 9^ 
is Borel-Leroy summable of order k in the distributional sense to F{g) for < g < R, with R > 0, 
if the following conditions are satisfied: 



(1) In accordance with Eq. |i76p , we set 

oo 



n=0 



(342) 



Then B(t), a priori holomorphic in some circle \t\ < A, admits a holomorphic continuation 
to the intersection of some neighbourhood of ]R,+ = {t G R : t > 0} with <C^ = {t € C : 
Imt > 0}. 

(a) The boundary value distributions i?(tibiO) = lim^^Q B^tzizie) exist\/t G ]R,+ and the following 
representation holds: 



1 poo , , 



l/k 



1/fc-l 



dt 



(343) 



for g belonging to the Nevanlinna disk of the g^^'' -plane, which is Cr = {5 G C : Re g > 
R-^], and where FP{Bj{t)) = ^ {B^ {t + iO) + B'^^\t + iO)). 



In this case, a criterion analogous to Theorem 1211 holds (see [189] for details). The generalization 
implied by Definition [22] is necessary, for instance, to prove the distributional Borel summability 
of the odd anharmonic oscillators x^"^"*"^ with m > 1. 



B.2 Distributional Borel Summability of the Stark Effect 

The Hydrogen Stark effect is observed when a single-electron atom is placed in an external constant 
electric field of strength F. If Z is the atomic number and /u the electron mass, then the Schrodinger 
equation becomes (with h = e = 1) 

■-^AiIj - —ip + Fxsj ip = Ftp, tp = 'ip{xi,X2,X3), r = Jxj + xl + xl, (344) 

where without loss of generality the electric field is directed along the X3 axis. Observe that in 
contrast to the discussion in Sec. IA.21 we are manifestly working in three dimensions here. The 
perturbation theory yields a subtle problem whose solution has direct connections to the odd 
anharmonic oscillators treated in Sec. 13. 5i In the presence of an electric field, the S0(4) symmetry 
of the hydrogen atom is broken; if however the field is constant the S0(2) symmetry around the 
direction of the field is conserved and the Schrodinger equation is separable in parabolic coordinates. 
Therefore the parabolic quantum numbers of the Hydrogen atom, ni, n2 and m, are used for 
the classification of the atomic states [237]. Explicitly, let us introduce the "squared" parabolic 
coordinates {u, v,(p), with < u,v < +00 and < (j) < 2-k: 

X2 

u = r + X3, V = r — X3, <f) = arctan — , (345a) 



inverted as 



xi = uvcos(p, X2 = uvsm(f>, X3 = -{u^ — v'^) . (345b) 
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Now, express the Laplace operator in the squared parabohc coordinate (see e.g. [237]); then exploit 
the cylindrical symmetry around the X3 axis to separate out the (p dependence writing f{u,v,(p) 
under the form f{u, v, cj)) = f{u, v)e^^'^ , m = 0, ±1, ±2, . . . and finally write f{u, v) under the form 
f{u,v) = ^{u)r]{v). Setting from now on = 1, we see that the Schrodinger equation p44p then 
reduces to a coupled system of ordinary differential equations: 

- l^Pr- + ""'"2^^^ ^(") - 2^^' + Fn^ ^{u) = Z,({u) (346) 

_ 1 d^ + ^(^) _ 2Eu^ r]{v) - F vS{v) = Z2r]{v) , (347) 

2 df^ 

where Zi + Z2 = Z and the boundary conditions ^(0) = 7?(0) = hold for ni = 0. The spectrum 
of (|344p is recovered in the following way: let Z^^'"^{E, F) with n,m € Nq be the eigenvalues of 
Eq. ([3l6]) . and Z^^''^{E,F) with n,m G Nq those of Eq. ([MT]) . Then the equation 

^nr,m^^^ F) + Z2"^''"(S, F)=Z (348) 

implicitly defines, for any triple (ni,n2,m), a family of functions E{ni,n2,'m, F) to be identified 
with the resonances of the Stark effect, which has no bound states as already recalled. In this context 
it should be remarked that the minus sign appearing in front of the quartic term in equation (j347p 
requires some care in the definition of the corresponding operator and thus of its eigenvalues. Note 
moreover that for F = and E < 0, the hydrogenic eigenvalues are recovered because in this 
case Z"i'™(£;,0) = y/'^2E {ni + |m| + 1/2) and ''""(£;, 0) = \/^2l'(n2 + |m| + 1/2). Hence the 
condition Zi + Z2 = Z yields 

E{ni,n2,m) = - — ■ , , , . ' (^^9) 

2(ni + n2 + \m\ + ly 

Thus, m is the standard magnetic quantum number; ni and n2 are called parabolic quantum 
numbers, and the principal quantum number n has the expression n = ni + 71-2 + \m\ + 1 in terms 
of the parabolic quantum numbers. 

As already mentioned, no hydrogen eigenvalue persists no matter how small the intensity F of the 
applied electric field is. We have indeed: 

Theorem 23 For F ^ the Schrodinger operator H{F) whose domain and action are defined as 
follows 

D{H{F))=C^{n''); H{F)u= (^-j-A-^ + Fx3^ u, ueD{H{F)), (350) 

is essentially self-adjoint, i.e. it has one and only one self-adjoint extension. Its spectrum (iden- 
tified with the spectrum of the self-adjoint extension) extends from —00 to +00 and is absolutely 
continuous. 



Proof. The first proof goes back to Titchmarsh, and relies on the separability in squared parabolic 
coordinates. A more recent and direct proof can be found in [322]. 

The absolute continuity of the spectrum in particular excludes the occurrence of eigenvalues em- 
bedded in the continuum. All eigenvalues of the hydrogen atom are thus expected to turn into 
resonances, and this physical intuition can actually be proven. 

Theorem 24 Let E{ni,n2, m) = E{ni,n2,m; 0) = — Z^/ (ni + 71-2 + ?ti + 1)^ be any eigenvalue of 
the hydrogen atom, classified by the parabolic quantum numbers ni > 0, n2 > and by the magnetic 
quantum number m € Nq. Then for \F\ ^ it turns into a resonance E{ni,n2,m] F), namely: 
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1. Let < argF < vr. Then the operator family H{F) defined by i350\) has discrete spectrum. 

2. There is B{ni,n2, m) > such that the eigenvalue E{ni,n2,m) of the hydrogen atom is stable 
as an eigenvalue E{ni,n2,in; F) of H{F) for \F\ < B{ni,n2,m), VO < argF < vr. 

3. The eigenvalues E{ni,n2,m; F) are holomorphic functions of F in the complex sector 
Q{ni,n2,m) = {F EC : |-F| < B{ni,n2,m) ; < argF < vr}. They admit an analytic 
continuation up to any sector {F G C : |F| < B^{ni,n2,m) ; — 7r/2 + e < argF < 37r/2 — e}, 
Ve > 0, through the real axis. 

4- For F £ a, F ^ 0, any function E{ni,n2,m; F) is complex valued, with the properties 

Re E{ni,n2,m; F + iO) = Re E{ni,n2,m; F — iO), (351a) 
ImF(ni, ?i2, m; F + iO) = — Im F(ni, n2, m; F — i 0) . (351b) 

Any complex-valued function E{ni,n2,m; F), F £ ]R is a resonance of the Stark effect. More 
specifically, ReF(ni, 712, m; F) is the location of the resonance, and ImF(ni, n2, F) its 
width, i.e. its inverse life-time. 



Proof See [322]. 

We remark that the functions E{ni,n2,m; F) are resonances of the Stark effect in the strongest 
possible sense; namely they verify the following properties 

1. They are second-sheet poles of the resolvent. This means the following: Consider the resolvent 
operator R{F, E) = [H{F) — E]~^. Since the spectrum of H{F) is the whole of R, the function 
R^{F,E) = (^|i?(F, F)^) is a holomorphic function of E in the upper half-plane ImF > 
for any ijj £ L'^. Then R^(F, E) has a meromorphic continuation to ImF < 0. The poles of 
the continuation coincide with the functions E(ni,n2,m; F). 

2. They are the eigenvalues of the complex dilated operator uniquely associated with H{F) by 
the Balslev-Combes analytic dilation procedure (see e.g. [44], XIII. 10). More precisely: for 
6* G R let 

{U{e)i>){x) = e3^/V(e^/^x); {U(e)~^^l^){x) = e-3^/V(e"''/^x) = (C/(-0)^)(x) (352) 

be the unitary dilation operator in L'^{B?). Let H{F, 0) = 11(9) H{F) 17(9)-^ be the unitary 
image of H{F), expressed by 

-26* y -9 

H{F, 9) = -^A - -—- + Fe^xs . (353) 
2/i r 

Then the operator family H{F, 9) can be analyzed also for 9 complex, < Im < vr. It can be 
proven that for such values of 9, the operator H{F,9) has discrete spectrum. Its eigenvalues 
do not depend on 9 and coincide with the functions E{ni,n2,m] F). 

3. They are second-sheet poles of the scattering matrix S{E). This means that the scattering 
matrix for the Stark effect S{E, F) exists, is unitary and has the same analyticity properties 
of the resolvent. Namely, it can be analytically continued from the upper half-plane ImF > 
to a meromorphic function to the lower half-plane ImF < 0. The poles of the continuation 
coincide with the functions E{ni,n2,m; F). 
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Proof. Assertions 1, 2, 3 and 4 of the theorem have been proven in [322] through the separabihty 
of the problem in squared parabohc coordinates recalled above, and in [326] by a direct analysis 
of the operator H{F,6). The complete equivalence between the two methods is proven in [206]. 
Assertion 5 is proven in [327]. An important issue about the equivalence between the separated 
Eqs. (j346p and (j347p and the operator H{F) for complex values of F is the interpretation of (j347p 
as the analytic continuation of (I346P from F > to —F = e"^F. Remark 3 is proven in [322]. 

Another point completely clarified in this example is the meaning of divergent perturbation theory. 
For the Stark effect, the perturbative expansion of the energy eigenvalue E{ni,n2,m, F) reads (see 
Eq. (59) of [328]), 

oo 

E{m,n2,m,F) ~ J] E^J^l^F^, (354) 

N=0 

where F is the electric field strength. 

Theorem 25 Consider the perturbation series (35^ for the resonance of quantum numbers 
{ni,n-2]m): 

oo 

E{n^ ,n2,m,F)^Yl ^n%m • (355) 
N=0 

Then: 

1. If < argF < vr, < |F| < B{ni,n2,m), the series is Borel summahle to the eigenvalue 
E{ni,n2,m, F) of H[F) which exists for such complex values of the field; 

2. If F ^ R, < |F| < B{{ni,n2,m) the series is Borel summable in the distributional sense to 
the resonance E{ni,n2,m, F). 

Proof. Assertion 1 is proven in [322], assertion 2 in [190]. 

Finally, we remark that the eigenvalues E(ni,n2,m, F), < argF < ir can be continued across the 
real axis, and that their continuation to F real from above yields the resonances. Since both the 
Borel sum and the analytic continuation are unique, we can say that the Borel summability along 
the complex directions of the field F is already sufficient to determine the resonances existing for 
F real. However, the direct summation of the original, real series to the resonances requires the 
distributional Borel procedure. 

The proof of the existence of resonances and of the (ordinary) Borel summability of their divergent 
perturbation expansion has been extended to the case of the Stark effect of an arbitrary atomic 
system, namely a system of quantum particles interacting through two-body Coulomb potentials 
and under the action of a constant electric field (see Refs. [329,330]). 

C Perturbation Series in Quantum Chromodynamics 

We would like to summarize the results of the investigations in Refs. [331-334] in the light of 
the concepts described in the current review. Typically, in quantum chromodynamics (QCD), the 
problem is to estimate the nonperturbative contribution to physical quantities V which, on the 
perturbative level, can be expressed as 

oo 

P-J^A„a^(Q), (356) 

n=0 



105 



where the strong couphng a'^{Q) is a running couphng constant that depends the momentum 
scale Q. Infrared (IR) renormalons lead to singularities along the integration path of the Laplace- 
Borel integral which have been amply discussed in Sec. [3j In the light of a resurgent expansion like 
Eq. ()252p . which has been conjectured for QCD perturbation series (see Ref. [186,187]), it is natural 
to associate the smallest term of a perturbation series with a fundamental uncertainty that has to 
be ascribed to the perturbative expansion due to nonperturbative effects. This prescription is also 
consistent with the numerical data in Table [TTl where the one-instanton series for the splitting of 
energy levels of the double well finds a natural barrier at the level of two-instanton shifts. 

The leading nonalternating divergence of perturbative expansions in QCD is commonly referred to 
as the leading IR-renormalon pole. Suppose that we can reformulate, using the (5 function of QCD, 

(-^) ~ (-'^° (aI^)) ^ [J^) ^ (^) ' 

where A denotes the Landau-pole and C corrects for the renormalization scheme. Suppose, further- 
more, that no is the smallest term of the perturbation series, and that the renormalon poles are 
located at positive integer k as, a. function of the argument of the Borel transform, defined accord- 
ing to Eq. (I173p . We then have as an estimate for the fundamental uncertainty of the perturbative 
expansion, 

no-l 



n=0 



A2e-^\*'" . / k 



Ec^'^' (^) ■ im 



k=l 

The nonperturbatively corrected perturbation series thus has the structure 



V = T.^r.s- + Y.C''^[^) . (359) 



n=0 k=l 



In some cases the k = 1 IR-renormalon is absent and the leading IR-renormalon appears at k — 2 
only, as for example in the case of the photon two-point-function [335] . From the vast number of 
related field-theoretical investigations, which profit from variants of the Borel method, we restrict 
ourselves to indicating the articles and books corresponding to Refs. [59,177,182,188,336-339]. 
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